diff --git a/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml b/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml index b6e2b8448..defd25867 100644 --- a/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml @@ -125,6 +125,12 @@ imap_lo_l1c_pset: Logical_source: imap_lo_l1c_pset Logical_source_description: IMAP Mission IMAP-Lo Instrument Level-1C Data +imap_lo_l1c_quickmap: + <<: *instrument_base + Data_type: L1C_quickmap>Level-1C QuickMap (Histogram Sky Map) + Logical_source: imap_lo_l1c_quickmap + Logical_source_description: IMAP Mission IMAP-Lo Instrument Level-1C Sky Map + # Global attributes for different sensors and sky tiling types and durations imap_lo_l2_enamap: <<: *instrument_base diff --git a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml index e56880dbe..3a8806aa5 100644 --- a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml @@ -131,7 +131,7 @@ met: spin_angle: CATDESC: Spin angle in degrees FIELDNAM: Spin Angle - FILLVAL: -1.0e31 + FILLVAL: -1.0e+31 FORMAT: F8.3 UNITS: deg VALIDMAX: 360.0 diff --git a/imap_processing/cdf/config/imap_lo_l1c_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1c_variable_attrs.yaml index d8d14122e..5b6589c90 100644 --- a/imap_processing/cdf/config/imap_lo_l1c_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1c_variable_attrs.yaml @@ -11,7 +11,7 @@ default_attrs: &default default_float32_attrs: &default_float32 <<: *default - FILLVAL: -1.0e31 + FILLVAL: -1.0e+31 FORMAT: F12.6 VALIDMAX: 3.4028235e+38 VALIDMIN: -3.4028235e+38 @@ -317,4 +317,211 @@ label_vector_HAE: CATDESC: Cartesian components (x,y,z) FIELDNAM: Cartesian components FORMAT: A5 - VAR_TYPE: metadata \ No newline at end of file + VAR_TYPE: metadata + +# --------------------------------------------------------------------------- +# imap_lo_l1c_quickmap sky-map variables (dims: esa_level, ecl_lat, ecl_lon). +# NaN in the systematic-error maps is written as the default_float32 FILLVAL +# (-1.0e31); readers should treat that value as fill. +# --------------------------------------------------------------------------- +flux_sys_err: + <<: *default_float32 + CATDESC: Hydrogen flux systematic error from asymmetric G-factor bounds (geometric mean) + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen flux systematic error + LABLAXIS: flux sys err + UNITS: "1/(cm^2 s sr keV)" + +flux_sys_err_upper: + <<: *default_float32 + CATDESC: Upper Hydrogen flux systematic error from the lower G-factor bound + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen flux systematic error upper + LABLAXIS: flux sys err hi + UNITS: "1/(cm^2 s sr keV)" + +flux_sys_err_lower: + <<: *default_float32 + CATDESC: Lower Hydrogen flux systematic error from the upper G-factor bound + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen flux systematic error lower + LABLAXIS: flux sys err lo + UNITS: "1/(cm^2 s sr keV)" + +flux_var_total: + <<: *default_float32 + CATDESC: Total Hydrogen flux variance (statistical plus systematic) + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen total flux variance + LABLAXIS: flux var total + UNITS: "1/(cm^2 s sr keV)^2" + +cosalpha: + <<: *default_float32 + CATDESC: RAM-direction projection factor sin(pivot) * sin(spin-angle) + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: RAM projection factor cos alpha + LABLAXIS: cos alpha + UNITS: " " + +# imap_lo_l1c_quickmap coordinates +esa_level: + CATDESC: ESA energy level + FIELDNAM: ESA level + FORMAT: I1 + LABLAXIS: ESA level + UNITS: " " + VALIDMAX: 7 + VALIDMIN: 1 + VAR_TYPE: support_data + +ecl_lat: + CATDESC: Ecliptic latitude bin center (ECLIPJ2000) + FIELDNAM: Ecliptic latitude + FORMAT: F8.3 + LABLAXIS: ecl lat + UNITS: "degrees" + VALIDMAX: 90.0 + VALIDMIN: -90.0 + VAR_TYPE: support_data + +ecl_lon: + CATDESC: Ecliptic longitude bin center (ECLIPJ2000) + FIELDNAM: Ecliptic longitude + FORMAT: F8.3 + LABLAXIS: ecl lon + UNITS: "degrees" + VALIDMAX: 360.0 + VALIDMIN: 0.0 + VAR_TYPE: support_data + +# imap_lo_l1c_quickmap sky-map data variables +counts: + <<: *default_float32 + CATDESC: Hydrogen counts accumulated per sky-map bin + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen counts + LABLAXIS: counts + UNITS: " " + +exposure: + <<: *default_float32 + CATDESC: Exposure time accumulated per sky-map bin + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Exposure time + LABLAXIS: exposure + UNITS: "s" + +rate: + <<: *default_float32 + CATDESC: Hydrogen count rate (counts / exposure) + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen rate + LABLAXIS: rate + UNITS: "1/s" + +rate_var: + <<: *default_float32 + CATDESC: Hydrogen count rate variance + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen rate variance + LABLAXIS: rate var + UNITS: "1/s^2" + +flux: + <<: *default_float32 + CATDESC: Hydrogen flux (rate / (G-factor * energy)) + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen flux + LABLAXIS: flux + UNITS: "1/(cm^2 s sr keV)" + +flux_var: + <<: *default_float32 + CATDESC: Hydrogen flux statistical variance + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen flux variance + LABLAXIS: flux var + UNITS: "1/(cm^2 s sr keV)^2" + +background_rate: + <<: *default_float32 + CATDESC: Hydrogen background rate + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen background rate + LABLAXIS: bg rate + UNITS: "1/s" + +background_rate_var: + <<: *default_float32 + CATDESC: Hydrogen background rate variance + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen background rate variance + LABLAXIS: bg rate var + UNITS: "1/s^2" + +background_flux: + <<: *default_float32 + CATDESC: Hydrogen background flux + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen background flux + LABLAXIS: bg flux + UNITS: "1/(cm^2 s sr keV)" + +background_flux_var: + <<: *default_float32 + CATDESC: Hydrogen background flux variance + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Hydrogen background flux variance + LABLAXIS: bg flux var + UNITS: "1/(cm^2 s sr keV)^2" + +signal_to_noise: + <<: *default_float32 + CATDESC: Signal-to-noise ratio (Hydrogen rate / background rate) + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Signal to noise + LABLAXIS: S/N + UNITS: " " + +signal_to_noise_var: + <<: *default_float32 + CATDESC: Signal-to-noise ratio variance + DEPEND_0: esa_level + DEPEND_1: ecl_lat + DEPEND_2: ecl_lon + FIELDNAM: Signal to noise variance + LABLAXIS: S/N var + UNITS: " " \ No newline at end of file diff --git a/imap_processing/cli.py b/imap_processing/cli.py index 324a46378..b0c5419c0 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -1304,23 +1304,41 @@ def do_processing( elif self.data_level == "l1c": data_dict = {} - anc_dependencies: list = dependencies.get_file_paths( - source="lo", data_type="ancillary" - ) - science_files = dependencies.get_file_paths(source="lo", descriptor="de") - science_files += dependencies.get_file_paths( - source="lo", data_type="l1b", descriptor="goodtimes" - ) - science_files += dependencies.get_file_paths( - source="lo", data_type="l1b", descriptor="bgrates" - ) - science_files += dependencies.get_file_paths( - source="lo", data_type="l1b", descriptor="histrates" - ) - for file in science_files: - dataset = load_cdf(file) - data_dict[dataset.attrs["Logical_source"]] = dataset - datasets = lo_l1c.lo_l1c(data_dict, anc_dependencies) + if self.descriptor == "pset": + anc_dependencies: list = dependencies.get_file_paths( + source="lo", data_type="ancillary" + ) + science_files = dependencies.get_file_paths( + source="lo", descriptor="de" + ) + science_files += dependencies.get_file_paths( + source="lo", data_type="l1b", descriptor="goodtimes" + ) + science_files += dependencies.get_file_paths( + source="lo", data_type="l1b", descriptor="bgrates" + ) + science_files += dependencies.get_file_paths( + source="lo", data_type="l1b", descriptor="histrates" + ) + for file in science_files: + dataset = load_cdf(file) + data_dict[dataset.attrs["Logical_source"]] = dataset + datasets = lo_l1c.lo_l1c(data_dict, anc_dependencies) + + elif self.descriptor == "quickmap": + science_files = ( + dependencies.get_file_paths(source="lo", descriptor="de") + + dependencies.get_file_paths(source="lo", descriptor="nhk") + + dependencies.get_file_paths(source="lo", descriptor="histrates") + + dependencies.get_file_paths(source="lo", descriptor="goodtimes") + + dependencies.get_file_paths(source="lo", descriptor="bgrates") + ) + + for file in science_files: + dataset = load_cdf(file) + data_dict[dataset.attrs["Logical_source"]] = dataset + + datasets = lo_l1c.lo_l1c_quickmap(data_dict) elif self.data_level == "l2": data_dict = {} diff --git a/imap_processing/lo/ancillary_data/imap_lo_bg-rates-anti-ram-overrides_20250901_v001.csv b/imap_processing/lo/ancillary_data/imap_lo_bg-rates-anti-ram-overrides_20250901_v001.csv new file mode 100755 index 000000000..22c0e1539 --- /dev/null +++ b/imap_processing/lo/ancillary_data/imap_lo_bg-rates-anti-ram-overrides_20250901_v001.csv @@ -0,0 +1,5 @@ +year,doy,counts/s +2026,62,0.0014 +2026,64,0.0 +2026,65,0.0 +2026,91,0.03 diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py index aa0e6e4a6..77a148c67 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -64,8 +64,8 @@ class LoConstants: } # Default background-rate thresholds [counts/s] when no pivot range matches. - THRESHOLD_BG_RATE_RAM_DEFAULT: float = 0.0175 - THRESHOLD_BG_RATE_ANTI_RAM_DEFAULT: float = 0.00875 + THRESHOLD_BG_RATE_RAM_DEFAULT: float = 0.028 + THRESHOLD_BG_RATE_ANTI_RAM_DEFAULT: float = 0.014 # Maximum time gap [s] between consecutive histogram epochs before treating them as # separate intervals. @@ -76,6 +76,47 @@ class LoConstants: # cycles are covered at interval edges. GOODTIME_PADDING: float = 2.0 + N_COLAT_BINS: int = 30 + + # The following are indexed by ESA level (0-indexed, ESA level = index + 1) + ESA_ENERGY: ClassVar[list[float]] = [ + 0.016, + 0.030, + 0.056, + 0.106, + 0.200, + 0.405, + 0.787, + 1.821, + ] + GEO_FACTOR: ClassVar[list[float]] = [ + 7.0e-5, + 7.9e-5, + 9.7e-5, + 11.2e-5, + 14.0e-5, + 17.7e-5, + 22.5e-5, + 6.721e-5, + ] + GEO_FACTOR_ERR: ClassVar[list[float]] = [ + 4.9e-5, + 5.5e-5, + 6.8e-5, + 3.0e-5, + 4.5e-5, + 2.0e-5, + 1.4e-5, + 6.721e-5, + ] + + # GEO_FACTOR/GEO_FACTOR_ERR above are the raw, pre-recalibration + # values; the quickmap product multiplies them by GEO_FACTOR_SCALE, and derives the + # asymmetric upper/lower G-factor bounds using the two scale factors below. + GEO_FACTOR_SCALE: float = 0.63529412 + GEO_FACTOR_SCALE_UPPER: float = 1.57407407 + GEO_FACTOR_SCALE_LOWER: float = 0.36728395 + # Star-sensor spin-angle binning offset (fractional bin-index shift used when # computing sample centers), keyed by the IFB star-sync housekeeping state # (ifb_ctrl_star_sync). Flight software 4.8 enabled star sync ("EN"), diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 1a92a605f..6b7d4a4c2 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -5,21 +5,25 @@ from enum import Enum import numpy as np +import pandas as pd import xarray as xr from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.ena_maps.utils.corrections import ( add_spacecraft_position_and_velocity_to_pset, ) +from imap_processing.lo.constants import LoConstants as c # noqa: N813 from imap_processing.spice.geometry import ( SpiceFrame, frame_transform_az_el, + get_spacecraft_to_instrument_spin_phase_offset, ) from imap_processing.spice.repoint import get_pointing_times_from_id from imap_processing.spice.spin import get_spin_number from imap_processing.spice.time import ( met_to_ttj2000ns, ttj2000ns_to_et, + ttj2000ns_to_met, ) N_ESA_ENERGY_STEPS = 7 @@ -36,6 +40,10 @@ OFF_ANGLE_BIN_EDGES = np.linspace(-2, 2, N_OFF_ANGLE_BINS + 1) OFF_ANGLE_BIN_CENTERS = (OFF_ANGLE_BIN_EDGES[:-1] + OFF_ANGLE_BIN_EDGES[1:]) / 2 +# Fill value written in place of NaN for quickmap float maps. Matches the +# default_float32 FILLVAL in imap_lo_l1c_variable_attrs.yaml. +QUICKMAP_FLOAT_FILLVAL = -1.0e31 + class FilterType(str, Enum): """ @@ -832,9 +840,10 @@ def compute_pointing_directions( dps_az_el = np.stack([spin, off], axis=-1) # Transform from DPS Az/El to the destination frame's lon/lat - return frame_transform_az_el( + az_el = frame_transform_az_el( et, dps_az_el, SpiceFrame.IMAP_DPS, to_frame, degrees=True ) + return np.asarray(az_el).reshape(len(spin_angles), len(off_angles), 2) def set_pointing_directions( @@ -877,3 +886,317 @@ def set_pointing_directions( dims=["epoch", "spin_angle", "off_angle"], attrs=attr_mgr.get_variable_attributes("hae_latitude"), ) + + +def lo_l1c_quickmap( # noqa: PLR0912 + sci_dependencies: dict, +) -> list[xr.Dataset]: + """ + Build Lo L1C quickmap datasets from L1B science dependencies. + + Determines the ecliptic sky pointing for each spin-angle bin via SPICE + (IMAP_DPS -> ECLIPJ2000 at the good-time midpoint), accumulates histogram + counts and exposure times within good-time intervals, projects them onto an + ecliptic sky grid, and derives flux and signal-to-noise maps for each ESA + energy level. + + Parameters + ---------- + sci_dependencies : dict + Dictionary of pre-computed L1B datasets keyed by product name. Must + contain ``"imap_lo_l1b_goodtimes"``, ``"imap_lo_l1b_bgrates"``, and + ``"imap_lo_l1b_histrates"``. + + Returns + ------- + list[xr.Dataset] + List of xarray Datasets, one per ESA energy level, each containing the + sky-map variables (counts, exposure, rates, flux, signal-to-noise, etc.) + on the ecliptic grid. + """ + attr_mgr = ImapCdfAttributes() + attr_mgr.add_instrument_global_attrs(instrument="lo") + attr_mgr.add_instrument_variable_attrs(instrument="lo", level="l1c") + + # Extract good-times and background rates from pre-computed dependencies + goodtimes_ds = sci_dependencies["imap_lo_l1b_goodtimes"] + bgrates_ds = sci_dependencies["imap_lo_l1b_bgrates"] + + pivot_angle = goodtimes_ds["pivot"].item() + gt_begin = goodtimes_ds["gt_start_met"].values + gt_end = goodtimes_ds["gt_end_met"].values + + # L1B histogram spin bins (spin_bin_6 <- L1A azimuth_6, which is just the + # onboard bin indices 0..59) are hardware spin-phase bins referenced to the + # spacecraft spin pulse -- NOT the instrument/DPS spin angle. Convert a bin + # center to the IMAP_DPS azimuth by ADDING the spacecraft->instrument spin-phase + # offset, exactly as the L1B star-sensor product does + # (calculate_star_sensor_profiles_by_group: + # ``spin_angle = start_angle_offset + sample_centers``). The L1B DE product + # needs no such offset because its spin_bin comes from a SPICE HAE->IMAP_DPS + # transform of the actual look direction, already in the instrument frame. + bin_width_deg = 360.0 / c.N_SPIN_ANGLE_BINS + sc_bin_centers = (np.arange(c.N_SPIN_ANGLE_BINS) + 0.5) * bin_width_deg + dps_offset_deg = ( + get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) * 360.0 + ) + dps_bin_angles = np.mod(sc_bin_centers + dps_offset_deg, 360.0) + + # Ecliptic sky pointing for each spin-angle bin, from SPICE at the good-time + # midpoint: pass the DPS azimuth (instrument spin angle) with a single boresight + # off-angle (0); the IMAP_DPS frame supplies the attitude. + # + # Validated: with the real IMAP_DPS/attitude CK furnished, this reproduces the + # legacy quickmap counts cell-for-cell (<= 2-count differences in 3 of 7 ESA + # levels, from sub-cell spin-axis rounding), confirming both the spin-phase + # offset and the IMAP_DPS azimuth zero-point match the legacy NEP convention. + # See tools/verify_lo_quickmap_equivalence.py. + pointing_epoch = met_to_ttj2000ns((gt_begin.min() + gt_end.max()) / 2.0) + az_el = compute_pointing_directions( + pointing_epoch, + pivot_angle, + spin_angles=dps_bin_angles, + off_angles=np.array([0.0]), + to_frame=SpiceFrame.ECLIPJ2000, + ) + ecl_lons = az_el[:, 0, 0] + ecl_lats = az_el[:, 0, 1] + + pivot_df = pd.DataFrame( + { + "bin_index": np.arange(c.N_SPIN_ANGLE_BINS), + "bins": dps_bin_angles, + "bin_ecl_lon": ecl_lons, + "bin_ecl_lat": ecl_lats, + } + ) + + # Filter histogram records to good-time intervals and accumulate per spin-angle bin + hist_ds = sci_dependencies["imap_lo_l1b_histrates"] + hist_met = ttj2000ns_to_met(hist_ds["epoch"].values) + mask = np.any( + (hist_met[:, np.newaxis] >= gt_begin) & (hist_met[:, np.newaxis] <= gt_end), + axis=1, + ) + + map_dataframes = [] + for esa_level in range(c.N_ESA_LEVELS): + hist_counts = np.sum(hist_ds["h_counts"].values[mask, esa_level, :].T, axis=1) + exposure = np.sum( + hist_ds["exposure_time_6deg"].values[mask, esa_level, :].T, axis=1 + ) + + esa_df = pd.DataFrame( + { + "bin_index": np.arange(c.N_SPIN_ANGLE_BINS), + "counts": hist_counts, + "expo": exposure, + } + ) + df = pivot_df.merge(esa_df, on="bin_index") + df.insert(0, "esa_level", esa_level + 1) + map_dataframes.append(df) + + map_df = pd.concat(map_dataframes, ignore_index=True) + + geo_scaled = [g * c.GEO_FACTOR_SCALE for g in c.GEO_FACTOR] + dg_scaled = [d * c.GEO_FACTOR_SCALE for d in c.GEO_FACTOR_ERR] + geo_err_upper = [ + float(np.hypot(g * (c.GEO_FACTOR_SCALE_UPPER - 1.0), d)) + for g, d in zip(geo_scaled, dg_scaled, strict=True) + ] + geo_err_lower = [ + float(np.hypot(g * (1.0 - c.GEO_FACTOR_SCALE_LOWER), d)) + for g, d in zip(geo_scaled, dg_scaled, strict=True) + ] + + # cosalpha uses the measured pivot with the +4 deg empirical offset + pivot_rad = np.radians(pivot_angle + 4.0) + + # Project onto N_COLAT_BINS x N_SPIN_ANGLE_BINS ecliptic sky grid and + # apply flux calibration + shape = (c.N_ESA_LEVELS, c.N_COLAT_BINS, c.N_SPIN_ANGLE_BINS) + ( + h_cnts_map, + exposure_map, + h_rate_map, + h_rate_var, + h_flux_map, + h_fvar_map, + h_fser_map, + h_fseu_map, + h_fsel_map, + h_fvto_map, + back_rate_map, + back_rate_var, + back_flux_map, + back_flux_var, + stonoise_map, + stonoise_var_map, + cosalpha_map, + ) = [np.zeros(shape) for _ in range(17)] + + for esa in range(c.N_ESA_LEVELS): + df = map_df[map_df["esa_level"] == esa + 1] + ps_ra = df["bin_ecl_lon"].values + ps_dec = df["bin_ecl_lat"].values + counts = df["counts"].values + expo_vals = df["expo"].values + + # Per-ESA calibration constants + energy = c.ESA_ENERGY[esa] + geo = geo_scaled[esa] + dgeu = geo_err_upper[esa] + dgel = geo_err_lower[esa] + # bgrates is shaped (epoch, esa_step); select the ESA step and take the + # single epoch to get a scalar (matches the pset path's `.values[0]`). + h_bgrate = float( + bgrates_ds["h_background_rates"].sel(esa_step=esa + 1).values[0] + ) + # The systematic flux error requires a positive lower G-factor bound + # (geo - dgel); otherwise it is undefined and written as FILLVAL. + sys_err_valid = geo > dgel + + for ia in range(c.N_SPIN_ANGLE_BINS): + imap = int(ps_ra[ia] * c.N_SPIN_ANGLE_BINS / 360.0) + if imap == c.N_SPIN_ANGLE_BINS: + imap = 0 + jmap = int((90.0 + ps_dec[ia]) * c.N_COLAT_BINS / 180.0) + if jmap == c.N_COLAT_BINS: + jmap = 0 + h_cnts_map[esa, jmap, imap] += counts[ia] + exposure_map[esa, jmap, imap] += expo_vals[ia] + + # RAM-direction projection factor sin(pivot) * sin(spin-angle), where + # the spin-angle is the instrument (DPS) bin angle carried in "bins". + alpha = np.radians(df["bins"].values[ia]) + cosalpha_map[esa, jmap, imap] = np.sin(pivot_rad) * np.sin(alpha) + + for imap in range(c.N_SPIN_ANGLE_BINS): + for jmap in range(c.N_COLAT_BINS): + expo = exposure_map[esa, jmap, imap] + if expo > 0: + back_rate_map[esa, jmap, imap] = h_bgrate + back_rate_var[esa, jmap, imap] = h_bgrate / expo + h_rate_map[esa, jmap, imap] = h_cnts_map[esa, jmap, imap] / expo + h_flux_map[esa, jmap, imap] = h_rate_map[esa, jmap, imap] / ( + geo * energy + ) + + # Asymmetric systematic flux error from the recalibrated + # G-factor bounds. The upper/lower flux excursions come from the + # lower/upper G-factor bounds respectively; h_fser is their + # geometric mean. Written as FILLVAL when the lower bound would + # drive the G-factor non-positive (systematic error undefined). + h_mid = h_flux_map[esa, jmap, imap] + if sys_err_valid: + h_hi = h_mid * geo / (geo - dgel) + h_lo = h_mid * geo / (geo + dgeu) + dh_hi = h_hi - h_mid + dh_lo = h_mid - h_lo + h_fser_map[esa, jmap, imap] = np.sqrt(dh_hi * dh_lo) + h_fseu_map[esa, jmap, imap] = dh_hi + h_fsel_map[esa, jmap, imap] = dh_lo + else: + h_fser_map[esa, jmap, imap] = QUICKMAP_FLOAT_FILLVAL + h_fseu_map[esa, jmap, imap] = QUICKMAP_FLOAT_FILLVAL + h_fsel_map[esa, jmap, imap] = QUICKMAP_FLOAT_FILLVAL + + back_flux_map[esa, jmap, imap] = h_bgrate / (geo * energy) + back_flux_var[esa, jmap, imap] = ( + back_rate_var[esa, jmap, imap] / (geo * energy) ** 2 + ) + h_rate_var[esa, jmap, imap] = h_rate_map[esa, jmap, imap] / expo + + if h_bgrate > 0.0: + stonoise_map[esa, jmap, imap] = ( + h_rate_map[esa, jmap, imap] / h_bgrate + ) + if ( + back_rate_map[esa, jmap, imap] > 0.0 + and h_rate_map[esa, jmap, imap] > 0.0 + ): + stonoise_var_map[esa, jmap, imap] = ( + h_rate_var[esa, jmap, imap] + / back_rate_map[esa, jmap, imap] ** 2 + + back_rate_var[esa, jmap, imap] + / h_rate_map[esa, jmap, imap] ** 2 + ) + + if h_cnts_map[esa, jmap, imap] > 0.0: + h_fvar_map[esa, jmap, imap] = ( + h_flux_map[esa, jmap, imap] ** 2 + / h_cnts_map[esa, jmap, imap] + ) + + # Total flux variance depends on the systematic error, so it is + # FILLVAL wherever the systematic error is undefined. + if not sys_err_valid: + h_fvto_map[esa, jmap, imap] = QUICKMAP_FLOAT_FILLVAL + elif h_cnts_map[esa, jmap, imap] > 0.0: + h_fvto_map[esa, jmap, imap] = ( + h_fvar_map[esa, jmap, imap] + + h_fser_map[esa, jmap, imap] ** 2 + ) + else: + h_fvto_map[esa, jmap, imap] = h_fser_map[esa, jmap, imap] ** 2 + + step_lon = 360.0 / c.N_SPIN_ANGLE_BINS + step_lat = 180.0 / c.N_COLAT_BINS + longitude_centers = np.arange(step_lon / 2, 360.0, step_lon) + ecl_lat_centers = np.arange(-90.0 + step_lat / 2, 90.0, step_lat) + + dims = ["esa_level", "ecl_lat", "ecl_lon"] + map_variables = { + "counts": h_cnts_map, + "exposure": exposure_map, + "rate": h_rate_map, + "rate_var": h_rate_var, + "flux": h_flux_map, + "flux_var": h_fvar_map, + "flux_sys_err": h_fser_map, + "flux_sys_err_upper": h_fseu_map, + "flux_sys_err_lower": h_fsel_map, + "flux_var_total": h_fvto_map, + "background_rate": back_rate_map, + "background_rate_var": back_rate_var, + "background_flux": back_flux_map, + "background_flux_var": back_flux_var, + "signal_to_noise": stonoise_map, + "signal_to_noise_var": stonoise_var_map, + "cosalpha": cosalpha_map, + } + data_vars = { + name: xr.DataArray( + data, + dims=dims, + attrs=attr_mgr.get_variable_attributes(name, check_schema=False), + ) + for name, data in map_variables.items() + } + + coords = { + "esa_level": xr.DataArray( + np.arange(1, c.N_ESA_LEVELS + 1), + dims=["esa_level"], + attrs=attr_mgr.get_variable_attributes("esa_level", check_schema=False), + ), + "ecl_lat": xr.DataArray( + ecl_lat_centers, + dims=["ecl_lat"], + attrs=attr_mgr.get_variable_attributes("ecl_lat", check_schema=False), + ), + "ecl_lon": xr.DataArray( + longitude_centers, + dims=["ecl_lon"], + attrs=attr_mgr.get_variable_attributes("ecl_lon", check_schema=False), + ), + } + + return [ + xr.Dataset( + data_vars, + coords=coords, + attrs=attr_mgr.get_global_attributes("imap_lo_l1c_quickmap"), + ) + ] diff --git a/imap_processing/tests/lo/test_lo_l1c.py b/imap_processing/tests/lo/test_lo_l1c.py index 40f6d4c08..b205b4624 100644 --- a/imap_processing/tests/lo/test_lo_l1c.py +++ b/imap_processing/tests/lo/test_lo_l1c.py @@ -6,18 +6,21 @@ from imap_processing import imap_module_directory from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes +from imap_processing.lo.constants import LoConstants from imap_processing.lo.l1c.lo_l1c import ( N_ESA_ENERGY_STEPS, N_OFF_ANGLE_BINS, N_SPIN_ANGLE_BINS, OFF_ANGLE_BIN_CENTERS, PSET_SHAPE, + QUICKMAP_FLOAT_FILLVAL, FilterType, calculate_exposure_times, compute_pointing_directions, create_pset_counts, filter_goodtimes, lo_l1c, + lo_l1c_quickmap, set_background_rates, set_pointing_directions, ) @@ -733,6 +736,235 @@ def test_set_pointing_directions_pivot_angle(attr_mgr, pivot_angle): np.testing.assert_allclose(actual_off_angles, expected_off_angles, rtol=1e-10) +# --------------------------------------------------------------------------- +# Integration test for the IMAP-Lo quickmap product (lo_l1c_quickmap). +# +# This exercises the whole quickmap path end-to-end on small synthetic L1B +# dependencies (good-times, bgrates, histrates) and asserts the structural and +# physical invariants the product must preserve. The sky pointing is obtained +# from SPICE (frame_transform_az_el), which is mocked here since a CK/attitude +# kernel is not available in the test environment. +# --------------------------------------------------------------------------- + +# Good-time window (MET seconds) that all "in-window" inputs fall inside. +_QM_GT_START = 511_000_000.0 +_QM_GT_END = 511_000_600.0 +_QM_PIVOT = 90.0 +# histrates epochs: first three MET values are inside the good-time window, +# the last two are outside it (and carry large counts/exposure that must be +# excluded by good-time filtering). +_QM_IN_METS = [511_000_150.0, 511_000_200.0, 511_000_250.0] +_QM_OUT_METS = [510_990_000.0, 511_010_000.0] + + +@pytest.fixture +def quickmap_inputs(): + """Small synthetic sci_dependencies for lo_l1c_quickmap.""" + n_esa = LoConstants.N_ESA_LEVELS # 7 + n_spin = 60 # L1B histogram spin bins (6 deg) + + mets = np.array(_QM_IN_METS + _QM_OUT_METS) + in_idx = np.array([0, 1, 2]) + out_idx = np.array([3, 4]) + + # h_counts / exposure: shape (n_epoch, n_esa, n_spin). + # In-window epochs get modest, per-esa, per-bin values so many sky cells end + # up populated; out-of-window epochs get large values that must be excluded. + rng = np.random.default_rng(42) + h_counts = np.zeros((mets.size, n_esa, n_spin)) + exposure = np.zeros((mets.size, n_esa, n_spin)) + for i in in_idx: + h_counts[i] = rng.integers(0, 4, size=(n_esa, n_spin)).astype(float) + exposure[i] = 10.0 * (np.arange(n_esa)[:, None] + 1) # positive everywhere + for i in out_idx: + h_counts[i] = 999.0 + exposure[i] = 999.0 + + histrates = xr.Dataset( + { + "h_counts": (["epoch", "esa_step", "spin_bin"], h_counts), + "exposure_time_6deg": (["epoch", "esa_step", "spin_bin"], exposure), + }, + coords={"epoch": met_to_ttj2000ns(mets)}, + ) + + goodtimes = xr.Dataset( + { + "pivot": ("epoch", [_QM_PIVOT]), + "gt_start_met": ("epoch", [_QM_GT_START]), + "gt_end_met": ("epoch", [_QM_GT_END]), + }, + coords={"epoch": [0]}, + ) + + h_bgrate = np.array([0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07], dtype=np.float64) + bgrates = xr.Dataset( + {"h_background_rates": (["epoch", "esa_step"], h_bgrate[np.newaxis, :])}, + coords={"epoch": [0], "esa_step": np.arange(1, n_esa + 1)}, + ) + + sci_dependencies = { + "imap_lo_l1b_histrates": histrates, + "imap_lo_l1b_goodtimes": goodtimes, + "imap_lo_l1b_bgrates": bgrates, + } + + # Expected per-esa totals from the in-window epochs only (the sky projection + # conserves the per-esa sum). + expected_counts = h_counts[in_idx].sum(axis=(0, 2)) + expected_exposure = exposure[in_idx].sum(axis=(0, 2)) + + return { + "sci_dependencies": sci_dependencies, + "h_bgrate": h_bgrate, + "expected_counts": expected_counts, + "expected_exposure": expected_exposure, + } + + +@pytest.fixture +def quickmap_result(quickmap_inputs): + """Run lo_l1c_quickmap once and share the result across assertions. + + ``frame_transform_az_el`` is mocked with an identity DPS->ecliptic transform + (spin angle -> longitude, off angle -> latitude) so the pointing is + deterministic without a CK/attitude kernel. Like the real function, the mock + squeezes the singleton off-angle dimension (returns ``(n_spin, 2)``), so the + (n_spin, n_off, 2) reshape in compute_pointing_directions is exercised. + """ + + def _identity_squeezed(et, az_el, *args, **kwargs): + return np.asarray(az_el)[:, 0, :] + + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=_identity_squeezed, + ): + datasets = lo_l1c_quickmap(quickmap_inputs["sci_dependencies"]) + assert len(datasets) == 1 + return datasets[0], quickmap_inputs + + +_QM_MAP_VARS = [ + "counts", + "exposure", + "rate", + "rate_var", + "flux", + "flux_var", + "flux_sys_err", + "flux_sys_err_upper", + "flux_sys_err_lower", + "flux_var_total", + "background_rate", + "background_rate_var", + "background_flux", + "background_flux_var", + "signal_to_noise", + "signal_to_noise_var", + "cosalpha", +] + + +def test_quickmap_structure(quickmap_result): + """The quickmap dataset has the expected dims, coords, vars, and source.""" + ds, _ = quickmap_result + n_esa = LoConstants.N_ESA_LEVELS + n_colat = LoConstants.N_COLAT_BINS + n_lon = LoConstants.N_SPIN_ANGLE_BINS + + assert ds.sizes == {"esa_level": n_esa, "ecl_lat": n_colat, "ecl_lon": n_lon} + assert ds.attrs["Logical_source"] == "imap_lo_l1c_quickmap" + + np.testing.assert_array_equal(ds["esa_level"].values, np.arange(1, n_esa + 1)) + assert ds["ecl_lat"].size == n_colat + assert ds["ecl_lon"].size == n_lon + + for var in _QM_MAP_VARS: + assert var in ds.data_vars, f"missing quickmap variable {var}" + assert ds[var].dims == ("esa_level", "ecl_lat", "ecl_lon") + + +def test_quickmap_count_conservation(quickmap_result): + """Total counts per ESA equal the summed in-window histrates counts.""" + ds, inp = quickmap_result + per_esa_total = ds["counts"].values.sum(axis=(1, 2)) + np.testing.assert_allclose(per_esa_total, inp["expected_counts"]) + + +def test_quickmap_exposure_conservation(quickmap_result): + """Total exposure per ESA equals the summed in-window histrates exposure.""" + ds, inp = quickmap_result + per_esa_total = ds["exposure"].values.sum(axis=(1, 2)) + np.testing.assert_allclose(per_esa_total, inp["expected_exposure"], rtol=1e-6) + + +def test_quickmap_goodtime_filtering(quickmap_result): + """Out-of-window epochs (counts=999/bin) are excluded by good-time filtering.""" + ds, inp = quickmap_result + # If the out-of-window epochs had leaked in, per-esa totals would be far + # larger than the in-window expectation (999 * 60 bins * 2 epochs). + per_esa_total = ds["counts"].values.sum(axis=(1, 2)) + np.testing.assert_allclose(per_esa_total, inp["expected_counts"]) + assert per_esa_total.max() < 999.0 * 60 + + +def test_quickmap_rate_and_flux_relationships(quickmap_result): + """Where exposed, rate = counts/exposure and flux = rate/(geo*energy).""" + ds, _ = quickmap_result + scale = LoConstants.GEO_FACTOR_SCALE + for e in range(LoConstants.N_ESA_LEVELS): + cnts = ds["counts"].values[e] + expo = ds["exposure"].values[e] + rate = ds["rate"].values[e] + flux = ds["flux"].values[e] + pos = expo > 0 + assert pos.any(), f"esa {e + 1} produced no exposed cells" + + np.testing.assert_allclose(rate[pos], cnts[pos] / expo[pos], rtol=1e-6) + assert np.all(rate[~pos] == 0.0) + + geo = LoConstants.GEO_FACTOR[e] * scale + energy = LoConstants.ESA_ENERGY[e] + np.testing.assert_allclose(flux[pos], rate[pos] / (geo * energy), rtol=1e-6) + + +def test_quickmap_background_rate(quickmap_result): + """Background rate map equals the input bgrate where exposed, else zero.""" + ds, inp = quickmap_result + for e in range(LoConstants.N_ESA_LEVELS): + expo = ds["exposure"].values[e] + brate = ds["background_rate"].values[e] + pos = expo > 0 + np.testing.assert_allclose(brate[pos], inp["h_bgrate"][e], rtol=1e-6) + assert np.all(brate[~pos] == 0.0) + + +def test_quickmap_systematic_error_fill(quickmap_result): + """Systematic-error fill follows the per-ESA G-factor validity condition. + + For each ESA the systematic flux error is only defined when the recalibrated + G-factor exceeds its lower error bound; otherwise the pipeline writes FILLVAL. + """ + ds, _ = quickmap_result + scale = LoConstants.GEO_FACTOR_SCALE + scale_lower = LoConstants.GEO_FACTOR_SCALE_LOWER + for e in range(LoConstants.N_ESA_LEVELS): + geo = LoConstants.GEO_FACTOR[e] * scale + dg = LoConstants.GEO_FACTOR_ERR[e] * scale + geo_err_lower = float(np.hypot(geo * (1.0 - scale_lower), dg)) + sys_err_valid = geo > geo_err_lower + + expo = ds["exposure"].values[e] + sys_err = ds["flux_sys_err"].values[e] + pos = expo > 0 + if sys_err_valid: + assert np.all(np.isfinite(sys_err[pos])) + assert np.all(sys_err[pos] != QUICKMAP_FLOAT_FILLVAL) + else: + assert np.all(sys_err[pos] == QUICKMAP_FLOAT_FILLVAL) + + def test_compute_pointing_directions_defaults(): """Default grid/frame reproduce the PSET (3600 x 40) IMAP_DPS->IMAP_HAE case.""" mock_et = 123456789.0