From 74b6aea7746d9a00d6ab7bb2078522d995fd8d21 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 17 Jul 2026 20:45:23 -0400 Subject: [PATCH 01/11] l1c quickmap implementation --- .../cdf/config/imap_lo_global_cdf_attrs.yaml | 6 + imap_processing/cli.py | 60 +++- imap_processing/lo/constants.py | 34 ++ imap_processing/lo/l1c/lo_l1c.py | 307 ++++++++++++++++++ 4 files changed, 390 insertions(+), 17 deletions(-) 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/cli.py b/imap_processing/cli.py index 6bacdde06..f14c35d5e 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -1300,23 +1300,49 @@ 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 + + quaternion_files = dependencies.get_file_paths( + source="spacecraft", descriptor="quaternions", data_type="l1a" + ) + quaternion_dependencies = [ + load_cdf(dep) for dep in list(set(quaternion_files)) + ] + quaternion_dependencies.sort(key=lambda ds: ds["epoch"].values[0]) + + datasets = lo_l1c.lo_l1c_quickmap(data_dict, quaternion_dependencies) elif self.data_level == "l2": data_dict = {} diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py index aa0e6e4a6..fa6c61fc7 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -76,6 +76,40 @@ 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, + ] + # 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 9d6b0161c..dc884efed 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -5,21 +5,29 @@ from enum import Enum import numpy as np +import pandas as pd import xarray as xr +from scipy.spatial.transform import Rotation 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.spacecraft.quaternions import assemble_quaternions from imap_processing.spice.geometry import ( SpiceFrame, + cartesian_to_spherical, frame_transform_az_el, + get_spacecraft_to_instrument_spin_phase_offset, + spherical_to_cartesian, ) 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 @@ -832,3 +840,302 @@ def set_pointing_directions( dims=["epoch", "spin_angle", "off_angle"], attrs=attr_mgr.get_variable_attributes("hae_latitude"), ) + + +def create_ra_dec( + spin_ecl_lon: float, spin_ecl_lat: float, pivot_angle: float +) -> tuple[list[float], list[float], list[float]]: + """ + Compute the ecliptic sky pointing for each of 60 spin-angle bins (6° each). + + All geometry is performed in ECLIPJ2000. For a spacecraft spin axis defined + by (spin_ecl_lon, spin_ecl_lat) and a pivot half-angle offset from that axis, + sweeps 360° in 60 equal steps and converts each pointing direction to ecliptic + longitude/latitude via ``cartesian_to_spherical``. The perpendicular plane is + oriented using the North Ecliptic Pole (= [0, 0, 1] in ECLIPJ2000) as a + reference, giving a frame with axes toward the NEP and toward the ram direction. + + Parameters + ---------- + spin_ecl_lon : float + Spin axis ecliptic longitude in degrees (ECLIPJ2000). + spin_ecl_lat : float + Spin axis ecliptic latitude in degrees (ECLIPJ2000). + pivot_angle : float + Half-angle offset from the spin axis in degrees. + + Returns + ------- + bin_centers : list of float + Spin angle bin centers (degrees) for each of the 60 bin centers. + ecl_lons : list of float + Ecliptic longitude (degrees, 0–360) for each of the 60 bin centers. + ecl_lats : list of float + Ecliptic latitude (degrees) for each of the 60 bin centers. + """ + pivot_angle_rad = np.radians(pivot_angle) + + # In ECLIPJ2000 the NEP is [0, 0, 1]. + nep_unit = np.array([0.0, 0.0, 1.0]) + spin_axis_unit = spherical_to_cartesian( + np.array([[1.0, spin_ecl_lon, spin_ecl_lat]]) + )[0] + spin_axis_unit /= np.linalg.norm(spin_axis_unit) + + # Build a right-handed frame in the plane perpendicular to the spin axis, + # anchored to the NEP so that spin-angle 0° points toward the pole. + spin_perp_toward_nep = nep_unit - np.dot(nep_unit, spin_axis_unit) * spin_axis_unit + spin_perp_toward_nep /= np.linalg.norm(spin_perp_toward_nep) + + spin_perp_toward_ram = np.cross(spin_axis_unit, spin_perp_toward_nep) + spin_perp_toward_ram /= np.linalg.norm(spin_perp_toward_ram) + + # Step through 60 spin-angle bin centers (3, 9, .., 357) and compute + # the 3D pointing direction for each bin using the pivot-angle cone equation, + # then read off ecliptic lon/lat via cartesian_to_spherical (ECLIPJ2000). + bin_centers: list[float] = np.arange(3.0, 360.0, 6.0).tolist() # type: ignore + + ecl_lons: list[float] = [] + ecl_lats: list[float] = [] + for spin_angle_deg in bin_centers: + pointing = ( + np.cos(pivot_angle_rad) * spin_axis_unit + + np.sin(pivot_angle_rad) + * np.cos(np.radians(spin_angle_deg)) + * spin_perp_toward_nep + + np.sin(pivot_angle_rad) + * np.sin(np.radians(spin_angle_deg)) + * spin_perp_toward_ram + ) + _, ecl_lon, ecl_lat = cartesian_to_spherical(pointing[None])[0] + ecl_lons.append(ecl_lon) + ecl_lats.append(ecl_lat) + + return bin_centers, ecl_lons, ecl_lats + + +def lo_l1c_quickmap( + sci_dependencies: dict, + quaternion_datasets: list[xr.Dataset], +) -> list[xr.Dataset]: + """ + Build Lo L1C quickmap datasets from L1B science dependencies and quaternions. + + Computes the mean spin-axis direction from quaternion data, determines the + ecliptic sky pointing for each spin-angle bin, 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"``. + quaternion_datasets : list[xr.Dataset] + List of spacecraft attitude quaternion datasets to be concatenated and + used for spin-axis estimation. + + 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() + h_bgrate = float(bgrates_ds["h_background_rates"].values[0]) + gt_begin = goodtimes_ds["gt_start_met"].values + gt_end = goodtimes_ds["gt_end_met"].values + + # Compute mean spin-axis direction in ECLIPJ2000 from quaternion data + quaternion_ds = xr.concat(quaternion_datasets, dim="epoch") + attitude_ds = assemble_quaternions(quaternion_ds) + + attitude_met = attitude_ds["epoch"].values + attitude_mask = np.any( + (attitude_met[:, None] >= gt_begin) & (attitude_met[:, None] <= gt_end), axis=1 + ) + attitude_ds = attitude_ds.isel(epoch=attitude_mask) + + quaternion_array = np.column_stack( + [ + attitude_ds["quat_x"], + attitude_ds["quat_y"], + attitude_ds["quat_z"], + attitude_ds["quat_s"], + ] + ) + mean_spin_axis = ( + Rotation.from_quat(quaternion_array).apply([0.0, 0.0, 1.0]).mean(axis=0) + ) + mean_spin_axis /= np.linalg.norm(mean_spin_axis) + spin_ecl_lon, spin_ecl_lat = cartesian_to_spherical(mean_spin_axis[None])[0, 1:] + spin_ecl_lon = float(spin_ecl_lon) + spin_ecl_lat = float(spin_ecl_lat) + + # Compute ecliptic sky pointing for each of 60 spin-angle bins + bins, ecl_lons, ecl_lats = create_ra_dec(spin_ecl_lon, spin_ecl_lat, pivot_angle) + pivot_df = pd.DataFrame( + {"bins": bins, "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[:, None] >= gt_begin) & (hist_met[:, None] <= gt_end), axis=1 + ) + + nep_roll = round( + get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) + * c.N_SPIN_ANGLE_BINS + ) + + 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 + ) + nep_counts = np.roll(hist_counts, nep_roll) + nep_exposure = np.roll(exposure, nep_roll) + + esa_df = pd.DataFrame( + {"bins": bins, "counts": nep_counts, "expo": nep_exposure} + ) + df = pivot_df.merge(esa_df, on="bins") + df.insert(0, "esa_level", esa_level + 1) + df["spin_ecl_lon"] = spin_ecl_lon + df["spin_ecl_lat"] = spin_ecl_lat + map_dataframes.append(df) + + map_df = pd.concat(map_dataframes, ignore_index=True) + + # 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_fvto_map, + back_rate_map, + back_rate_var, + back_flux_map, + back_flux_var, + stonoise_map, + stonoise_var_map, + ) = [np.zeros(shape) for _ in range(14)] + + 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 + + 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] + + 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: + energy = c.ESA_ENERGY[esa] + geo = c.GEO_FACTOR[esa] + dge = c.GEO_FACTOR_ERR[esa] + + 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 + ) + h_fser_map[esa, jmap, imap] = ( + h_rate_map[esa, jmap, imap] * dge / (geo**2 * energy) + ) + 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] + ) + 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"] + return [ + xr.Dataset( + { + "counts": xr.DataArray(h_cnts_map, dims=dims), + "exposure": xr.DataArray(exposure_map, dims=dims), + "rate": xr.DataArray(h_rate_map, dims=dims), + "rate_var": xr.DataArray(h_rate_var, dims=dims), + "flux": xr.DataArray(h_flux_map, dims=dims), + "flux_var": xr.DataArray(h_fvar_map, dims=dims), + "flux_sys_err": xr.DataArray(h_fser_map, dims=dims), + "flux_var_total": xr.DataArray(h_fvto_map, dims=dims), + "background_rate": xr.DataArray(back_rate_map, dims=dims), + "background_rate_var": xr.DataArray(back_rate_var, dims=dims), + "background_flux": xr.DataArray(back_flux_map, dims=dims), + "background_flux_var": xr.DataArray(back_flux_var, dims=dims), + "signal_to_noise": xr.DataArray(stonoise_map, dims=dims), + "signal_to_noise_var": xr.DataArray(stonoise_var_map, dims=dims), + }, + coords={ + "esa_level": np.arange(1, c.N_ESA_LEVELS + 1), + "ecl_lat": ecl_lat_centers, + "ecl_lon": longitude_centers, + }, + attrs=attr_mgr.get_global_attributes("imap_lo_l1c_quickmap"), + ) + ] From 419601829d2a6612f770dcbf398d353b63884424 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 10:20:17 -0400 Subject: [PATCH 02/11] minor tweaks to get the code in line with map_SCFrame.py --- .../imap_lo_bg-rates-anti-ram-overrides_20250901_v001.csv | 5 +++++ imap_processing/lo/constants.py | 4 ++-- imap_processing/lo/l1c/lo_l1c.py | 4 +++- 3 files changed, 10 insertions(+), 3 deletions(-) create mode 100755 imap_processing/lo/ancillary_data/imap_lo_bg-rates-anti-ram-overrides_20250901_v001.csv 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 fa6c61fc7..7d4ea8295 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. diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index dc884efed..855c0ddf4 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -952,7 +952,6 @@ def lo_l1c_quickmap( bgrates_ds = sci_dependencies["imap_lo_l1b_bgrates"] pivot_angle = goodtimes_ds["pivot"].item() - h_bgrate = float(bgrates_ds["h_background_rates"].values[0]) gt_begin = goodtimes_ds["gt_start_met"].values gt_end = goodtimes_ds["gt_end_met"].values @@ -1064,6 +1063,9 @@ def lo_l1c_quickmap( energy = c.ESA_ENERGY[esa] geo = c.GEO_FACTOR[esa] dge = c.GEO_FACTOR_ERR[esa] + h_bgrate = float( + bgrates_ds["h_background_rates"].sel(esa_step=esa + 1) + ) back_rate_map[esa, jmap, imap] = h_bgrate back_rate_var[esa, jmap, imap] = h_bgrate / expo From bc6ebfa2ce3884eef56df85ed0fa67d7bf3e4eed Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 14:43:26 -0400 Subject: [PATCH 03/11] quickmap product after validation --- .../config/imap_lo_l1b_variable_attrs.yaml | 2 +- .../config/imap_lo_l1c_variable_attrs.yaml | 211 +++++++++++++++++- imap_processing/lo/constants.py | 7 + imap_processing/lo/l1c/lo_l1c.py | 153 ++++++++++--- 4 files changed, 337 insertions(+), 36 deletions(-) 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..4cb398fe4 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: " " + VALIDMIN: 1 + VALIDMAX: 7 + VAR_TYPE: support_data + +ecl_lat: + CATDESC: Ecliptic latitude bin center (ECLIPJ2000) + FIELDNAM: Ecliptic latitude + FORMAT: F8.3 + LABLAXIS: ecl lat + UNITS: "degrees" + VALIDMIN: -90.0 + VALIDMAX: 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" + VALIDMIN: 0.0 + VALIDMAX: 360.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/lo/constants.py b/imap_processing/lo/constants.py index 7d4ea8295..77a148c67 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -110,6 +110,13 @@ class LoConstants: 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 855c0ddf4..10bd8b4fb 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -44,6 +44,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): """ @@ -914,7 +918,7 @@ def create_ra_dec( return bin_centers, ecl_lons, ecl_lats -def lo_l1c_quickmap( +def lo_l1c_quickmap( # noqa: PLR0912 sci_dependencies: dict, quaternion_datasets: list[xr.Dataset], ) -> list[xr.Dataset]: @@ -1019,6 +1023,24 @@ def lo_l1c_quickmap( map_df = pd.concat(map_dataframes, ignore_index=True) + # June 6, 2026 geometric-factor recalibration (see map_SCFrame_V2). The raw + # G-factors are rescaled, and asymmetric upper/lower G-factor errors are formed + # by combining the multiplicative-bound offset with the raw per-step error in + # quadrature. + 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 (map_SCFrame_V2) + 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) @@ -1030,6 +1052,8 @@ def lo_l1c_quickmap( h_flux_map, h_fvar_map, h_fser_map, + h_fseu_map, + h_fsel_map, h_fvto_map, back_rate_map, back_rate_var, @@ -1037,7 +1061,8 @@ def lo_l1c_quickmap( back_flux_var, stonoise_map, stonoise_var_map, - ) = [np.zeros(shape) for _ in range(14)] + 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] @@ -1046,6 +1071,20 @@ def lo_l1c_quickmap( 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: @@ -1056,26 +1095,41 @@ def lo_l1c_quickmap( 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 NEP-frame bin center already 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: - energy = c.ESA_ENERGY[esa] - geo = c.GEO_FACTOR[esa] - dge = c.GEO_FACTOR_ERR[esa] - h_bgrate = float( - bgrates_ds["h_background_rates"].sel(esa_step=esa + 1) - ) - 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 ) - h_fser_map[esa, jmap, imap] = ( - h_rate_map[esa, jmap, imap] * dge / (geo**2 * 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 @@ -1102,6 +1156,12 @@ def lo_l1c_quickmap( 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 @@ -1115,29 +1175,56 @@ def lo_l1c_quickmap( 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( - { - "counts": xr.DataArray(h_cnts_map, dims=dims), - "exposure": xr.DataArray(exposure_map, dims=dims), - "rate": xr.DataArray(h_rate_map, dims=dims), - "rate_var": xr.DataArray(h_rate_var, dims=dims), - "flux": xr.DataArray(h_flux_map, dims=dims), - "flux_var": xr.DataArray(h_fvar_map, dims=dims), - "flux_sys_err": xr.DataArray(h_fser_map, dims=dims), - "flux_var_total": xr.DataArray(h_fvto_map, dims=dims), - "background_rate": xr.DataArray(back_rate_map, dims=dims), - "background_rate_var": xr.DataArray(back_rate_var, dims=dims), - "background_flux": xr.DataArray(back_flux_map, dims=dims), - "background_flux_var": xr.DataArray(back_flux_var, dims=dims), - "signal_to_noise": xr.DataArray(stonoise_map, dims=dims), - "signal_to_noise_var": xr.DataArray(stonoise_var_map, dims=dims), - }, - coords={ - "esa_level": np.arange(1, c.N_ESA_LEVELS + 1), - "ecl_lat": ecl_lat_centers, - "ecl_lon": longitude_centers, - }, + data_vars, + coords=coords, attrs=attr_mgr.get_global_attributes("imap_lo_l1c_quickmap"), ) ] From b8efcb3fcf12b05aba0cf05d386d75a58cfffa57 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 14:52:29 -0400 Subject: [PATCH 04/11] integration tests before refactor --- imap_processing/tests/lo/test_lo_l1c.py | 239 ++++++++++++++++++++++++ 1 file changed, 239 insertions(+) diff --git a/imap_processing/tests/lo/test_lo_l1c.py b/imap_processing/tests/lo/test_lo_l1c.py index 01c1f46c3..15fb65dee 100644 --- a/imap_processing/tests/lo/test_lo_l1c.py +++ b/imap_processing/tests/lo/test_lo_l1c.py @@ -6,17 +6,20 @@ 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, create_pset_counts, filter_goodtimes, lo_l1c, + lo_l1c_quickmap, set_background_rates, set_pointing_directions, ) @@ -729,3 +732,239 @@ def test_set_pointing_directions_pivot_angle(attr_mgr, pivot_angle): # dps_az_el[:, :, 1] should have the adjusted off angles repeated across spin actual_off_angles = dps_az_el[0, :, 1] # Take first spin 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 + quaternions and asserts the structural and physical invariants +# the product must preserve. It is intended as a behaviour lock-in *before* +# refactoring lo_l1c_quickmap: a correct refactor must keep all of these true. +# --------------------------------------------------------------------------- + +# 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] + + +def _make_quaternion_ds(): + """Build a raw 10 Hz quaternion dataset as assemble_quaternions expects. + + Every 10 Hz sample carries the same fixed quaternion (a 45-degree rotation + about x), so the mean spin axis is well defined and not aligned with the + ecliptic pole (which would make create_ra_dec degenerate). The packet times + are inside the good-time window so the attitude mask is non-empty. + """ + times = np.arange(511_000_100.0, 511_000_120.0, 1.0) # 20 packets, in-window + quat = [np.sin(np.radians(22.5)), 0.0, 0.0, np.cos(np.radians(22.5))] + data_vars = {"sciencedata1hz_quat_10_hz_time": ("epoch", times)} + for quat_i, comp in enumerate(quat): + for i in range(10): + data_vars[f"fsw_acs_quat_10_hz_buffered_{i + quat_i * 10}"] = ( + "epoch", + np.full(times.size, comp), + ) + return xr.Dataset(data_vars, coords={"epoch": times}) + + +@pytest.fixture +def quickmap_inputs(): + """Small synthetic sci_dependencies + quaternions 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 (roll + projection + # both conserve 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, + "quaternions": [_make_quaternion_ds()], + "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.""" + datasets = lo_l1c_quickmap( + quickmap_inputs["sci_dependencies"], quickmap_inputs["quaternions"] + ) + 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) From 6000261c0910dc16746c0a292a7498c4c4fa8f40 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 15:20:44 -0400 Subject: [PATCH 05/11] computing each histogram bin's sky pointing at its true (bin_center + offset) angle, instead of a np.roll after the fact --- imap_processing/lo/l1c/lo_l1c.py | 90 +++++++++++++++++++++----------- 1 file changed, 59 insertions(+), 31 deletions(-) diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 10bd8b4fb..fa780edb3 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -847,14 +847,17 @@ def set_pointing_directions( def create_ra_dec( - spin_ecl_lon: float, spin_ecl_lat: float, pivot_angle: float + spin_ecl_lon: float, + spin_ecl_lat: float, + pivot_angle: float, + spin_angles: np.ndarray | None = None, ) -> tuple[list[float], list[float], list[float]]: """ - Compute the ecliptic sky pointing for each of 60 spin-angle bins (6° each). + Compute the ecliptic sky pointing for a set of spin-angle bins. All geometry is performed in ECLIPJ2000. For a spacecraft spin axis defined by (spin_ecl_lon, spin_ecl_lat) and a pivot half-angle offset from that axis, - sweeps 360° in 60 equal steps and converts each pointing direction to ecliptic + each spin-angle is converted to a pointing direction and then to ecliptic longitude/latitude via ``cartesian_to_spherical``. The perpendicular plane is oriented using the North Ecliptic Pole (= [0, 0, 1] in ECLIPJ2000) as a reference, giving a frame with axes toward the NEP and toward the ram direction. @@ -867,15 +870,18 @@ def create_ra_dec( Spin axis ecliptic latitude in degrees (ECLIPJ2000). pivot_angle : float Half-angle offset from the spin axis in degrees. + spin_angles : numpy.ndarray, optional + Spin-angle bin centers (degrees) in the NEP-anchored frame to evaluate. + Defaults to the 60 centers of the 6° histogram bins (3, 9, .., 357). Returns ------- bin_centers : list of float - Spin angle bin centers (degrees) for each of the 60 bin centers. + The spin-angle bin centers (degrees) that were evaluated. ecl_lons : list of float - Ecliptic longitude (degrees, 0–360) for each of the 60 bin centers. + Ecliptic longitude (degrees, 0–360) for each bin center. ecl_lats : list of float - Ecliptic latitude (degrees) for each of the 60 bin centers. + Ecliptic latitude (degrees) for each bin center. """ pivot_angle_rad = np.radians(pivot_angle) @@ -894,10 +900,12 @@ def create_ra_dec( spin_perp_toward_ram = np.cross(spin_axis_unit, spin_perp_toward_nep) spin_perp_toward_ram /= np.linalg.norm(spin_perp_toward_ram) - # Step through 60 spin-angle bin centers (3, 9, .., 357) and compute - # the 3D pointing direction for each bin using the pivot-angle cone equation, - # then read off ecliptic lon/lat via cartesian_to_spherical (ECLIPJ2000). - bin_centers: list[float] = np.arange(3.0, 360.0, 6.0).tolist() # type: ignore + # For each spin-angle bin center compute the 3D pointing direction using the + # pivot-angle cone equation, then read off ecliptic lon/lat via + # cartesian_to_spherical (ECLIPJ2000). + if spin_angles is None: + spin_angles = np.arange(3.0, 360.0, 6.0) + bin_centers: list[float] = np.asarray(spin_angles, dtype=float).tolist() ecl_lons: list[float] = [] ecl_lats: list[float] = [] @@ -911,7 +919,7 @@ def create_ra_dec( * np.sin(np.radians(spin_angle_deg)) * spin_perp_toward_ram ) - _, ecl_lon, ecl_lat = cartesian_to_spherical(pointing[None])[0] + _, ecl_lon, ecl_lat = cartesian_to_spherical(pointing[np.newaxis])[0] ecl_lons.append(ecl_lon) ecl_lats.append(ecl_lat) @@ -965,7 +973,9 @@ def lo_l1c_quickmap( # noqa: PLR0912 attitude_met = attitude_ds["epoch"].values attitude_mask = np.any( - (attitude_met[:, None] >= gt_begin) & (attitude_met[:, None] <= gt_end), axis=1 + (attitude_met[:, np.newaxis] >= gt_begin) + & (attitude_met[:, np.newaxis] <= gt_end), + axis=1, ) attitude_ds = attitude_ds.isel(epoch=attitude_mask) @@ -981,26 +991,46 @@ def lo_l1c_quickmap( # noqa: PLR0912 Rotation.from_quat(quaternion_array).apply([0.0, 0.0, 1.0]).mean(axis=0) ) mean_spin_axis /= np.linalg.norm(mean_spin_axis) - spin_ecl_lon, spin_ecl_lat = cartesian_to_spherical(mean_spin_axis[None])[0, 1:] + spin_ecl_lon, spin_ecl_lat = cartesian_to_spherical(mean_spin_axis[np.newaxis])[ + 0, 1: + ] spin_ecl_lon = float(spin_ecl_lon) spin_ecl_lat = float(spin_ecl_lat) - # Compute ecliptic sky pointing for each of 60 spin-angle bins - bins, ecl_lons, ecl_lats = create_ra_dec(spin_ecl_lon, spin_ecl_lat, pivot_angle) + # Rotate the spacecraft-frame histogram spin-bin centers into the NEP-anchored + # frame that `create_ra_dec` expects. The spacecraft->instrument spin-phase + # offset is applied continuously as an angle. + 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 + nep_offset_deg = ( + get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) * 360.0 + ) + # The offset is ADDED (positive) to the spacecraft spin-bin angle to get the + # NEP-frame angle. Sign convention comes from + # get_spacecraft_to_instrument_spin_phase_offset (angle measured positively from + # the S/C +x-axis, per the imap_130 frames kernel) combined with the right-handed + # NEP->RAM sweep in create_ra_dec. + nep_bin_angles = np.mod(sc_bin_centers + nep_offset_deg, 360.0) + + # Compute ecliptic sky pointing for each spin-angle bin at its NEP-frame angle + bins, ecl_lons, ecl_lats = create_ra_dec( + spin_ecl_lon, spin_ecl_lat, pivot_angle, spin_angles=nep_bin_angles + ) pivot_df = pd.DataFrame( - {"bins": bins, "bin_ecl_lon": ecl_lons, "bin_ecl_lat": ecl_lats} + { + "bin_index": np.arange(c.N_SPIN_ANGLE_BINS), + "bins": bins, + "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[:, None] >= gt_begin) & (hist_met[:, None] <= gt_end), axis=1 - ) - - nep_roll = round( - get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) - * c.N_SPIN_ANGLE_BINS + (hist_met[:, np.newaxis] >= gt_begin) & (hist_met[:, np.newaxis] <= gt_end), + axis=1, ) map_dataframes = [] @@ -1009,13 +1039,15 @@ def lo_l1c_quickmap( # noqa: PLR0912 exposure = np.sum( hist_ds["exposure_time_6deg"].values[mask, esa_level, :].T, axis=1 ) - nep_counts = np.roll(hist_counts, nep_roll) - nep_exposure = np.roll(exposure, nep_roll) esa_df = pd.DataFrame( - {"bins": bins, "counts": nep_counts, "expo": nep_exposure} + { + "bin_index": np.arange(c.N_SPIN_ANGLE_BINS), + "counts": hist_counts, + "expo": exposure, + } ) - df = pivot_df.merge(esa_df, on="bins") + df = pivot_df.merge(esa_df, on="bin_index") df.insert(0, "esa_level", esa_level + 1) df["spin_ecl_lon"] = spin_ecl_lon df["spin_ecl_lat"] = spin_ecl_lat @@ -1023,10 +1055,6 @@ def lo_l1c_quickmap( # noqa: PLR0912 map_df = pd.concat(map_dataframes, ignore_index=True) - # June 6, 2026 geometric-factor recalibration (see map_SCFrame_V2). The raw - # G-factors are rescaled, and asymmetric upper/lower G-factor errors are formed - # by combining the multiplicative-bound offset with the raw per-step error in - # quadrature. 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 = [ @@ -1038,7 +1066,7 @@ def lo_l1c_quickmap( # noqa: PLR0912 for g, d in zip(geo_scaled, dg_scaled, strict=True) ] - # cosalpha uses the measured pivot with the +4 deg empirical offset (map_SCFrame_V2) + # 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 From 4d635d4744c59fa40a19973175e99bae4b5bfaba Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 16:05:30 -0400 Subject: [PATCH 06/11] Extract reusable compute_pointing_directions from set_pointing_directions Generalize the Lo L1C pointing-direction geometry so it can be reused by other products (e.g. the l1c quickmap) without duplicating the DPS az/el -> sky transform. - Add compute_pointing_directions(epoch, pivot_angle, spin_angles=None, off_angles=None, to_frame=IMAP_HAE): the pure-geometry core returning the raw (n_spin, n_off, 2) lon/lat array for an arbitrary spin/off-angle grid and destination frame. - set_pointing_directions now delegates to it with the PSET defaults (3600x40 grid, IMAP_DPS -> IMAP_HAE) and keeps its DataArray output, so the PSET path is unchanged. - Add unit tests for the defaults, custom grid/frame, and the delegation. Co-Authored-By: Claude Opus 4.8 (1M context) --- imap_processing/lo/l1c/lo_l1c.py | 71 ++++++++++++++++++---- imap_processing/tests/lo/test_lo_l1c.py | 79 +++++++++++++++++++++++++ 2 files changed, 137 insertions(+), 13 deletions(-) diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 9d6b0161c..1a92a605f 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -780,6 +780,63 @@ def set_background_rates( return bg_rates_data, bg_stat_uncert_data, bg_sys_err_data +def compute_pointing_directions( + epoch: float, + pivot_angle: float, + spin_angles: np.ndarray | None = None, + off_angles: np.ndarray | None = None, + to_frame: SpiceFrame = SpiceFrame.IMAP_HAE, +) -> np.ndarray: + """ + Transform a DPS spin/off-angle grid to sky longitude/latitude at an epoch. + + For the given ``epoch`` the instrument despun (IMAP_DPS) spin-angle / off-angle + grid is transformed into ``to_frame`` longitude/latitude using SPICE. Off-angles + are measured relative to ``pivot_angle`` (elevation = 90 - pivot_angle + off). + + This is the reusable geometry core shared by pointing-direction products: pass + the desired spin/off-angle bin centers and destination frame. + + Parameters + ---------- + epoch : float + The epoch time in TTJ2000ns. + pivot_angle : float + The pivot angle in degrees. + Off-angles are adjusted relative to this pivot angle before transformation. + spin_angles : numpy.ndarray, optional + Spin-angle bin centers in degrees. Defaults to ``SPIN_ANGLE_BIN_CENTERS`` + (the 3600-bin PSET grid). + off_angles : numpy.ndarray, optional + Off-angle bin centers in degrees. Defaults to ``OFF_ANGLE_BIN_CENTERS`` + (the 40-bin PSET grid). + to_frame : SpiceFrame, optional + Destination reference frame. Defaults to ``SpiceFrame.IMAP_HAE``. + + Returns + ------- + numpy.ndarray + Array of shape ``(n_spin, n_off, 2)`` where ``[..., 0]`` is longitude and + ``[..., 1]`` is latitude, both in degrees in ``to_frame``. + """ + if spin_angles is None: + spin_angles = SPIN_ANGLE_BIN_CENTERS + if off_angles is None: + off_angles = OFF_ANGLE_BIN_CENTERS + + et = ttj2000ns_to_et(epoch) + # create a meshgrid of spin and off angles using the bin centers + spin, off = np.meshgrid(spin_angles, off_angles, indexing="ij") + # off_angles need to account for the pivot_angle + off = off + (90 - pivot_angle) + 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( + et, dps_az_el, SpiceFrame.IMAP_DPS, to_frame, degrees=True + ) + + def set_pointing_directions( epoch: float, attr_mgr: ImapCdfAttributes, @@ -809,19 +866,7 @@ def set_pointing_directions( hae_latitude : xr.DataArray The HAE latitude for each spin and off angle bin. """ - et = ttj2000ns_to_et(epoch) - # create a meshgrid of spin and off angles using the bin centers - spin, off = np.meshgrid( - SPIN_ANGLE_BIN_CENTERS, OFF_ANGLE_BIN_CENTERS, indexing="ij" - ) - # off_angles need to account for the pivot_angle - off += 90 - pivot_angle - dps_az_el = np.stack([spin, off], axis=-1) - - # Transform from DPS Az/El to HAE lon/lat - hae_az_el = frame_transform_az_el( - et, dps_az_el, SpiceFrame.IMAP_DPS, SpiceFrame.IMAP_HAE, degrees=True - ) + hae_az_el = compute_pointing_directions(epoch, pivot_angle) return xr.DataArray( data=hae_az_el[np.newaxis, :, :, 0].astype(np.float64), diff --git a/imap_processing/tests/lo/test_lo_l1c.py b/imap_processing/tests/lo/test_lo_l1c.py index 01c1f46c3..40f6d4c08 100644 --- a/imap_processing/tests/lo/test_lo_l1c.py +++ b/imap_processing/tests/lo/test_lo_l1c.py @@ -14,12 +14,14 @@ PSET_SHAPE, FilterType, calculate_exposure_times, + compute_pointing_directions, create_pset_counts, filter_goodtimes, lo_l1c, set_background_rates, set_pointing_directions, ) +from imap_processing.spice.geometry import SpiceFrame from imap_processing.spice.time import met_to_ttj2000ns @@ -729,3 +731,80 @@ def test_set_pointing_directions_pivot_angle(attr_mgr, pivot_angle): # dps_az_el[:, :, 1] should have the adjusted off angles repeated across spin actual_off_angles = dps_az_el[0, :, 1] # Take first spin angle np.testing.assert_allclose(actual_off_angles, expected_off_angles, rtol=1e-10) + + +def test_compute_pointing_directions_defaults(): + """Default grid/frame reproduce the PSET (3600 x 40) IMAP_DPS->IMAP_HAE case.""" + mock_et = 123456789.0 + mock_az_el = np.stack( + np.meshgrid(np.arange(3600), np.arange(40), indexing="ij"), axis=-1 + ) + with ( + patch("imap_processing.lo.l1c.lo_l1c.ttj2000ns_to_et") as mock_ttj2000ns_to_et, + patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el" + ) as mock_frame_transform, + ): + mock_ttj2000ns_to_et.return_value = mock_et + mock_frame_transform.return_value = mock_az_el + + result = compute_pointing_directions(1000000000.0, 90) + + # Returns the raw (n_spin, n_off, 2) array, not a DataArray. + assert result.shape == (3600, 40, 2) + call_args = mock_frame_transform.call_args + assert call_args[0][1].shape == (3600, 40, 2) # dps_az_el grid + assert call_args[0][2] == SpiceFrame.IMAP_DPS # from_frame + assert call_args[0][3] == SpiceFrame.IMAP_HAE # default to_frame + + +def test_compute_pointing_directions_custom_grid_and_frame(): + """Custom spin/off angles and destination frame are honored.""" + spin_angles = np.arange(3.0, 360.0, 6.0) # 60 bins + off_angles = np.array([0.0]) # single pivot-cone off-angle + pivot_angle = 75.0 + with ( + patch("imap_processing.lo.l1c.lo_l1c.ttj2000ns_to_et") as mock_ttj2000ns_to_et, + patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el" + ) as mock_frame_transform, + ): + mock_ttj2000ns_to_et.return_value = 123456789.0 + mock_frame_transform.side_effect = lambda et, az_el, *a, **k: az_el + + result = compute_pointing_directions( + 1000000000.0, + pivot_angle, + spin_angles=spin_angles, + off_angles=off_angles, + to_frame=SpiceFrame.ECLIPJ2000, + ) + + assert result.shape == (60, 1, 2) + # Spin component matches the requested spin angles. + np.testing.assert_allclose(result[:, 0, 0], spin_angles) + # Off component is the single off-angle offset by (90 - pivot_angle). + np.testing.assert_allclose(result[:, 0, 1], 90 - pivot_angle) + # Destination frame is forwarded. + assert mock_frame_transform.call_args[0][3] == SpiceFrame.ECLIPJ2000 + + +def test_set_pointing_directions_delegates(attr_mgr): + """set_pointing_directions wraps compute_pointing_directions output unchanged.""" + mock_az_el = np.stack( + np.meshgrid(np.arange(3600), np.arange(40), indexing="ij"), axis=-1 + ).astype(float) + with patch( + "imap_processing.lo.l1c.lo_l1c.compute_pointing_directions" + ) as mock_compute: + mock_compute.return_value = mock_az_el + + hae_longitude, hae_latitude = set_pointing_directions( + 1000000000.0, attr_mgr, 90 + ) + + mock_compute.assert_called_once_with(1000000000.0, 90) + assert hae_longitude.dims == ("epoch", "spin_angle", "off_angle") + assert hae_longitude.shape == (1, 3600, 40) + np.testing.assert_array_equal(hae_longitude.values[0], mock_az_el[:, :, 0]) + np.testing.assert_array_equal(hae_latitude.values[0], mock_az_el[:, :, 1]) From 1c0ab585d568103d5ecb58ac3c9cbc3867ee7d32 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 16:30:01 -0400 Subject: [PATCH 07/11] lo_l1c_quickmap: get sky pointing from SPICE via compute_pointing_directions Replace the quaternion-derived mean-spin-axis + analytic create_ra_dec pointing with a SPICE IMAP_DPS -> ECLIPJ2000 transform at the good-time midpoint, using the reusable compute_pointing_directions helper. - Drop the quaternion_datasets parameter; lo_l1c_quickmap now takes only sci_dependencies. Update the CLI quickmap branch accordingly (no longer loads spacecraft quaternion files). - Remove create_ra_dec and the assemble_quaternions / Rotation / cartesian_to_spherical / spherical_to_cartesian imports it needed. - Keep nep_bin_angles solely for the cosalpha (RAM projection) factor. - Integration tests: drop the synthetic quaternion dataset and mock frame_transform_az_el with an identity DPS->ecliptic transform, since a CK/attitude kernel is not available in the test environment. Caveat documented in-code: routing through IMAP_DPS assumes the L1B histogram spin-bin index is the DPS spin angle (mounting handled by the frame kernel). This has not been validated against a known-good sky map; if the histogram is in raw spacecraft spin phase the map would be rotated by the spin-phase offset. Co-Authored-By: Claude Opus 4.8 (1M context) --- imap_processing/cli.py | 10 +- imap_processing/lo/l1c/lo_l1c.py | 170 +++++------------------- imap_processing/tests/lo/test_lo_l1c.py | 49 +++---- 3 files changed, 52 insertions(+), 177 deletions(-) diff --git a/imap_processing/cli.py b/imap_processing/cli.py index f14c35d5e..680c7ab46 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -1334,15 +1334,7 @@ def do_processing( dataset = load_cdf(file) data_dict[dataset.attrs["Logical_source"]] = dataset - quaternion_files = dependencies.get_file_paths( - source="spacecraft", descriptor="quaternions", data_type="l1a" - ) - quaternion_dependencies = [ - load_cdf(dep) for dep in list(set(quaternion_files)) - ] - quaternion_dependencies.sort(key=lambda ds: ds["epoch"].values[0]) - - datasets = lo_l1c.lo_l1c_quickmap(data_dict, quaternion_dependencies) + datasets = lo_l1c.lo_l1c_quickmap(data_dict) elif self.data_level == "l2": data_dict = {} diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index dab2fbc4f..56f5c3510 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -7,20 +7,16 @@ import numpy as np import pandas as pd import xarray as xr -from scipy.spatial.transform import Rotation 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.spacecraft.quaternions import assemble_quaternions from imap_processing.spice.geometry import ( SpiceFrame, - cartesian_to_spherical, frame_transform_az_el, get_spacecraft_to_instrument_spin_phase_offset, - spherical_to_cartesian, ) from imap_processing.spice.repoint import get_pointing_times_from_id from imap_processing.spice.spin import get_spin_number @@ -891,97 +887,17 @@ def set_pointing_directions( ) -def create_ra_dec( - spin_ecl_lon: float, - spin_ecl_lat: float, - pivot_angle: float, - spin_angles: np.ndarray | None = None, -) -> tuple[list[float], list[float], list[float]]: - """ - Compute the ecliptic sky pointing for a set of spin-angle bins. - - All geometry is performed in ECLIPJ2000. For a spacecraft spin axis defined - by (spin_ecl_lon, spin_ecl_lat) and a pivot half-angle offset from that axis, - each spin-angle is converted to a pointing direction and then to ecliptic - longitude/latitude via ``cartesian_to_spherical``. The perpendicular plane is - oriented using the North Ecliptic Pole (= [0, 0, 1] in ECLIPJ2000) as a - reference, giving a frame with axes toward the NEP and toward the ram direction. - - Parameters - ---------- - spin_ecl_lon : float - Spin axis ecliptic longitude in degrees (ECLIPJ2000). - spin_ecl_lat : float - Spin axis ecliptic latitude in degrees (ECLIPJ2000). - pivot_angle : float - Half-angle offset from the spin axis in degrees. - spin_angles : numpy.ndarray, optional - Spin-angle bin centers (degrees) in the NEP-anchored frame to evaluate. - Defaults to the 60 centers of the 6° histogram bins (3, 9, .., 357). - - Returns - ------- - bin_centers : list of float - The spin-angle bin centers (degrees) that were evaluated. - ecl_lons : list of float - Ecliptic longitude (degrees, 0–360) for each bin center. - ecl_lats : list of float - Ecliptic latitude (degrees) for each bin center. - """ - pivot_angle_rad = np.radians(pivot_angle) - - # In ECLIPJ2000 the NEP is [0, 0, 1]. - nep_unit = np.array([0.0, 0.0, 1.0]) - spin_axis_unit = spherical_to_cartesian( - np.array([[1.0, spin_ecl_lon, spin_ecl_lat]]) - )[0] - spin_axis_unit /= np.linalg.norm(spin_axis_unit) - - # Build a right-handed frame in the plane perpendicular to the spin axis, - # anchored to the NEP so that spin-angle 0° points toward the pole. - spin_perp_toward_nep = nep_unit - np.dot(nep_unit, spin_axis_unit) * spin_axis_unit - spin_perp_toward_nep /= np.linalg.norm(spin_perp_toward_nep) - - spin_perp_toward_ram = np.cross(spin_axis_unit, spin_perp_toward_nep) - spin_perp_toward_ram /= np.linalg.norm(spin_perp_toward_ram) - - # For each spin-angle bin center compute the 3D pointing direction using the - # pivot-angle cone equation, then read off ecliptic lon/lat via - # cartesian_to_spherical (ECLIPJ2000). - if spin_angles is None: - spin_angles = np.arange(3.0, 360.0, 6.0) - bin_centers: list[float] = np.asarray(spin_angles, dtype=float).tolist() - - ecl_lons: list[float] = [] - ecl_lats: list[float] = [] - for spin_angle_deg in bin_centers: - pointing = ( - np.cos(pivot_angle_rad) * spin_axis_unit - + np.sin(pivot_angle_rad) - * np.cos(np.radians(spin_angle_deg)) - * spin_perp_toward_nep - + np.sin(pivot_angle_rad) - * np.sin(np.radians(spin_angle_deg)) - * spin_perp_toward_ram - ) - _, ecl_lon, ecl_lat = cartesian_to_spherical(pointing[np.newaxis])[0] - ecl_lons.append(ecl_lon) - ecl_lats.append(ecl_lat) - - return bin_centers, ecl_lons, ecl_lats - - def lo_l1c_quickmap( # noqa: PLR0912 sci_dependencies: dict, - quaternion_datasets: list[xr.Dataset], ) -> list[xr.Dataset]: """ - Build Lo L1C quickmap datasets from L1B science dependencies and quaternions. + Build Lo L1C quickmap datasets from L1B science dependencies. - Computes the mean spin-axis direction from quaternion data, determines the - ecliptic sky pointing for each spin-angle bin, 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. + 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 ---------- @@ -989,9 +905,6 @@ def lo_l1c_quickmap( # noqa: PLR0912 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"``. - quaternion_datasets : list[xr.Dataset] - List of spacecraft attitude quaternion datasets to be concatenated and - used for spin-axis estimation. Returns ------- @@ -1012,59 +925,44 @@ def lo_l1c_quickmap( # noqa: PLR0912 gt_begin = goodtimes_ds["gt_start_met"].values gt_end = goodtimes_ds["gt_end_met"].values - # Compute mean spin-axis direction in ECLIPJ2000 from quaternion data - quaternion_ds = xr.concat(quaternion_datasets, dim="epoch") - attitude_ds = assemble_quaternions(quaternion_ds) - - attitude_met = attitude_ds["epoch"].values - attitude_mask = np.any( - (attitude_met[:, np.newaxis] >= gt_begin) - & (attitude_met[:, np.newaxis] <= gt_end), - axis=1, - ) - attitude_ds = attitude_ds.isel(epoch=attitude_mask) - - quaternion_array = np.column_stack( - [ - attitude_ds["quat_x"], - attitude_ds["quat_y"], - attitude_ds["quat_z"], - attitude_ds["quat_s"], - ] - ) - mean_spin_axis = ( - Rotation.from_quat(quaternion_array).apply([0.0, 0.0, 1.0]).mean(axis=0) - ) - mean_spin_axis /= np.linalg.norm(mean_spin_axis) - spin_ecl_lon, spin_ecl_lat = cartesian_to_spherical(mean_spin_axis[np.newaxis])[ - 0, 1: - ] - spin_ecl_lon = float(spin_ecl_lon) - spin_ecl_lat = float(spin_ecl_lat) - - # Rotate the spacecraft-frame histogram spin-bin centers into the NEP-anchored - # frame that `create_ra_dec` expects. The spacecraft->instrument spin-phase - # offset is applied continuously as an angle. + # Histogram spin-bin centers (spacecraft/instrument spin phase, degrees). 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 + + # NEP-anchored spin angle, used only by the cosalpha (RAM projection) factor + # below. The spacecraft->instrument offset is ADDED (positive) to the spin-bin + # angle; sign convention comes from get_spacecraft_to_instrument_spin_phase_offset + # (angle measured positively from the S/C +x-axis, per the imap_130 frames kernel). nep_offset_deg = ( get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) * 360.0 ) - # The offset is ADDED (positive) to the spacecraft spin-bin angle to get the - # NEP-frame angle. Sign convention comes from - # get_spacecraft_to_instrument_spin_phase_offset (angle measured positively from - # the S/C +x-axis, per the imap_130 frames kernel) combined with the right-handed - # NEP->RAM sweep in create_ra_dec. nep_bin_angles = np.mod(sc_bin_centers + nep_offset_deg, 360.0) - # Compute ecliptic sky pointing for each spin-angle bin at its NEP-frame angle - bins, ecl_lons, ecl_lats = create_ra_dec( - spin_ecl_lon, spin_ecl_lat, pivot_angle, spin_angles=nep_bin_angles + # Ecliptic sky pointing for each spin-angle bin, from SPICE at the good-time + # midpoint. The instrument despun (IMAP_DPS) frame carries the spacecraft-> + # instrument mounting and attitude, so the raw spin-bin centers are passed as + # the DPS azimuth (spin angle) with a single boresight off-angle (0). + # + # NOTE (unvalidated): this assumes the L1B histogram spin-bin index is the + # IMAP_DPS spin angle (mounting handled by the frame kernel). It has not been + # checked against a known-good sky map -- if the histogram convention is raw + # spacecraft spin phase instead, the resulting map would be rotated by + # nep_offset_deg. + pointing_epoch = met_to_ttj2000ns((gt_begin.min() + gt_end.max()) / 2.0) + az_el = compute_pointing_directions( + pointing_epoch, + pivot_angle, + spin_angles=sc_bin_centers, + 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": bins, + "bins": nep_bin_angles, "bin_ecl_lon": ecl_lons, "bin_ecl_lat": ecl_lats, } @@ -1094,8 +992,6 @@ def lo_l1c_quickmap( # noqa: PLR0912 ) df = pivot_df.merge(esa_df, on="bin_index") df.insert(0, "esa_level", esa_level + 1) - df["spin_ecl_lon"] = spin_ecl_lon - df["spin_ecl_lat"] = spin_ecl_lat map_dataframes.append(df) map_df = pd.concat(map_dataframes, ignore_index=True) diff --git a/imap_processing/tests/lo/test_lo_l1c.py b/imap_processing/tests/lo/test_lo_l1c.py index 062f325d4..ecf3405b5 100644 --- a/imap_processing/tests/lo/test_lo_l1c.py +++ b/imap_processing/tests/lo/test_lo_l1c.py @@ -740,9 +740,10 @@ def test_set_pointing_directions_pivot_angle(attr_mgr, pivot_angle): # 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 + quaternions and asserts the structural and physical invariants -# the product must preserve. It is intended as a behaviour lock-in *before* -# refactoring lo_l1c_quickmap: a correct refactor must keep all of these true. +# 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. @@ -756,29 +757,9 @@ def test_set_pointing_directions_pivot_angle(attr_mgr, pivot_angle): _QM_OUT_METS = [510_990_000.0, 511_010_000.0] -def _make_quaternion_ds(): - """Build a raw 10 Hz quaternion dataset as assemble_quaternions expects. - - Every 10 Hz sample carries the same fixed quaternion (a 45-degree rotation - about x), so the mean spin axis is well defined and not aligned with the - ecliptic pole (which would make create_ra_dec degenerate). The packet times - are inside the good-time window so the attitude mask is non-empty. - """ - times = np.arange(511_000_100.0, 511_000_120.0, 1.0) # 20 packets, in-window - quat = [np.sin(np.radians(22.5)), 0.0, 0.0, np.cos(np.radians(22.5))] - data_vars = {"sciencedata1hz_quat_10_hz_time": ("epoch", times)} - for quat_i, comp in enumerate(quat): - for i in range(10): - data_vars[f"fsw_acs_quat_10_hz_buffered_{i + quat_i * 10}"] = ( - "epoch", - np.full(times.size, comp), - ) - return xr.Dataset(data_vars, coords={"epoch": times}) - - @pytest.fixture def quickmap_inputs(): - """Small synthetic sci_dependencies + quaternions for lo_l1c_quickmap.""" + """Small synthetic sci_dependencies for lo_l1c_quickmap.""" n_esa = LoConstants.N_ESA_LEVELS # 7 n_spin = 60 # L1B histogram spin bins (6 deg) @@ -828,14 +809,13 @@ def quickmap_inputs(): "imap_lo_l1b_bgrates": bgrates, } - # Expected per-esa totals from the in-window epochs only (roll + projection - # both conserve the per-esa sum). + # 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, - "quaternions": [_make_quaternion_ds()], "h_bgrate": h_bgrate, "expected_counts": expected_counts, "expected_exposure": expected_exposure, @@ -844,10 +824,17 @@ def quickmap_inputs(): @pytest.fixture def quickmap_result(quickmap_inputs): - """Run lo_l1c_quickmap once and share the result across assertions.""" - datasets = lo_l1c_quickmap( - quickmap_inputs["sci_dependencies"], quickmap_inputs["quaternions"] - ) + """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. + """ + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=lambda et, az_el, *args, **kwargs: az_el, + ): + datasets = lo_l1c_quickmap(quickmap_inputs["sci_dependencies"]) assert len(datasets) == 1 return datasets[0], quickmap_inputs From 9944e30ae8231edd70bfd19c32d5030e3979b16d Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 16:39:30 -0400 Subject: [PATCH 08/11] lo_l1c_quickmap: pass instrument (DPS) spin angle, resolving the offset convention The L1B histogram spin bins (spin_bin_6 <- L1A azimuth_6 = onboard bin indices 0..59) are hardware spin-phase bins referenced to the spacecraft spin pulse, not the instrument/DPS spin angle. The L1B star-sensor product converts the same onboard angle to the instrument frame by adding the spacecraft->instrument spin-phase offset (spin_angle = start_angle_offset + sample_centers); the L1B DE product needs no offset only because its spin_bin comes from a SPICE HAE->IMAP_DPS transform of the actual look direction. Accordingly, feed compute_pointing_directions the offset-applied DPS azimuth (sc_bin_centers + offset) rather than the raw bin centers, matching the legacy quickmap (which rolled the histogram by the same offset). Update the in-code note from "unvalidated assumption" to the resolved convention, leaving only a sky-map cross-check of the DPS azimuth zero-point as a follow-up. Co-Authored-By: Claude Opus 4.8 (1M context) --- imap_processing/lo/l1c/lo_l1c.py | 38 ++++++++++++++++---------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 56f5c3510..3234bf894 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -925,34 +925,34 @@ def lo_l1c_quickmap( # noqa: PLR0912 gt_begin = goodtimes_ds["gt_start_met"].values gt_end = goodtimes_ds["gt_end_met"].values - # Histogram spin-bin centers (spacecraft/instrument spin phase, degrees). + # 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 - - # NEP-anchored spin angle, used only by the cosalpha (RAM projection) factor - # below. The spacecraft->instrument offset is ADDED (positive) to the spin-bin - # angle; sign convention comes from get_spacecraft_to_instrument_spin_phase_offset - # (angle measured positively from the S/C +x-axis, per the imap_130 frames kernel). - nep_offset_deg = ( + dps_offset_deg = ( get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) * 360.0 ) - nep_bin_angles = np.mod(sc_bin_centers + nep_offset_deg, 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. The instrument despun (IMAP_DPS) frame carries the spacecraft-> - # instrument mounting and attitude, so the raw spin-bin centers are passed as - # the DPS azimuth (spin angle) with a single boresight off-angle (0). + # midpoint: pass the DPS azimuth (instrument spin angle) with a single boresight + # off-angle (0); the IMAP_DPS frame supplies the attitude. # - # NOTE (unvalidated): this assumes the L1B histogram spin-bin index is the - # IMAP_DPS spin angle (mounting handled by the frame kernel). It has not been - # checked against a known-good sky map -- if the histogram convention is raw - # spacecraft spin phase instead, the resulting map would be rotated by - # nep_offset_deg. + # NOTE: the spin-phase offset above is grounded in the L1B code, but a sky-map + # cross-check is still worthwhile to confirm the IMAP_DPS azimuth zero-point + # matches the legacy NEP-anchored convention beyond this offset. pointing_epoch = met_to_ttj2000ns((gt_begin.min() + gt_end.max()) / 2.0) az_el = compute_pointing_directions( pointing_epoch, pivot_angle, - spin_angles=sc_bin_centers, + spin_angles=dps_bin_angles, off_angles=np.array([0.0]), to_frame=SpiceFrame.ECLIPJ2000, ) @@ -962,7 +962,7 @@ def lo_l1c_quickmap( # noqa: PLR0912 pivot_df = pd.DataFrame( { "bin_index": np.arange(c.N_SPIN_ANGLE_BINS), - "bins": nep_bin_angles, + "bins": dps_bin_angles, "bin_ecl_lon": ecl_lons, "bin_ecl_lat": ecl_lats, } @@ -1065,7 +1065,7 @@ def lo_l1c_quickmap( # noqa: PLR0912 exposure_map[esa, jmap, imap] += expo_vals[ia] # RAM-direction projection factor sin(pivot) * sin(spin-angle), where - # the spin-angle is the NEP-frame bin center already carried in "bins". + # 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) From ab3f030785b28402e3b7803da0a1f8884b2fe592 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 16:58:22 -0400 Subject: [PATCH 09/11] lo_l1c_quickmap: fix pointing shape for singleton off-angle; note CK validation frame_transform_az_el squeezes singleton leading dims, so compute_pointing_directions returned (n_spin, 2) rather than (n_spin, 1, 2) when called with a single boresight off-angle (as the quickmap does). Reshape back to the documented (n_spin, n_off, 2) contract so az_el[:, 0, :] indexing works with real SPICE (the identity mock in the tests did not reproduce the squeeze). Also upgrade the in-code note: with the real IMAP_DPS/attitude CK furnished (tools/verify_lo_quickmap_equivalence.py), the SPICE pointing reproduces the legacy quickmap counts cell-for-cell (<= 2-count differences in 3 of 7 ESA levels, from sub-cell spin-axis rounding), confirming the spin-phase offset and the IMAP_DPS azimuth zero-point match the legacy NEP convention. Co-Authored-By: Claude Opus 4.8 (1M context) --- imap_processing/lo/l1c/lo_l1c.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 3234bf894..2b4bdab94 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -839,10 +839,13 @@ def compute_pointing_directions( off = off + (90 - pivot_angle) 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( + # Transform from DPS Az/El to the destination frame's lon/lat. + # frame_transform_az_el squeezes singleton leading dims, so reshape back to + # the documented (n_spin, n_off, 2) contract (e.g. a single off-angle bin). + 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( @@ -945,9 +948,11 @@ def lo_l1c_quickmap( # noqa: PLR0912 # midpoint: pass the DPS azimuth (instrument spin angle) with a single boresight # off-angle (0); the IMAP_DPS frame supplies the attitude. # - # NOTE: the spin-phase offset above is grounded in the L1B code, but a sky-map - # cross-check is still worthwhile to confirm the IMAP_DPS azimuth zero-point - # matches the legacy NEP-anchored convention beyond this offset. + # 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, From 0852d559f62c5f638b38820a6d688b32fabbbc64 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 20 Jul 2026 16:59:12 -0400 Subject: [PATCH 10/11] test: make quickmap pointing mock reproduce frame_transform_az_el squeeze The identity mock previously returned (n_spin, 1, 2), hiding the singleton-dim squeeze that real frame_transform_az_el performs. Return (n_spin, 2) instead so the reshape in compute_pointing_directions is actually exercised and a regression would fail the test. Co-Authored-By: Claude Opus 4.8 (1M context) --- imap_processing/tests/lo/test_lo_l1c.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/imap_processing/tests/lo/test_lo_l1c.py b/imap_processing/tests/lo/test_lo_l1c.py index ecf3405b5..b205b4624 100644 --- a/imap_processing/tests/lo/test_lo_l1c.py +++ b/imap_processing/tests/lo/test_lo_l1c.py @@ -828,11 +828,17 @@ def quickmap_result(quickmap_inputs): ``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. + 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=lambda et, az_el, *args, **kwargs: az_el, + side_effect=_identity_squeezed, ): datasets = lo_l1c_quickmap(quickmap_inputs["sci_dependencies"]) assert len(datasets) == 1 From f6abb154f80c17217d67804ab25304704965f352 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 24 Jul 2026 14:16:03 -0400 Subject: [PATCH 11/11] test fixes --- imap_processing/cdf/config/imap_lo_l1c_variable_attrs.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 4cb398fe4..5b6589c90 100644 --- a/imap_processing/cdf/config/imap_lo_l1c_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1c_variable_attrs.yaml @@ -381,8 +381,8 @@ esa_level: FORMAT: I1 LABLAXIS: ESA level UNITS: " " - VALIDMIN: 1 VALIDMAX: 7 + VALIDMIN: 1 VAR_TYPE: support_data ecl_lat: @@ -391,8 +391,8 @@ ecl_lat: FORMAT: F8.3 LABLAXIS: ecl lat UNITS: "degrees" - VALIDMIN: -90.0 VALIDMAX: 90.0 + VALIDMIN: -90.0 VAR_TYPE: support_data ecl_lon: @@ -401,8 +401,8 @@ ecl_lon: FORMAT: F8.3 LABLAXIS: ecl lon UNITS: "degrees" - VALIDMIN: 0.0 VALIDMAX: 360.0 + VALIDMIN: 0.0 VAR_TYPE: support_data # imap_lo_l1c_quickmap sky-map data variables