A placeholder quickmap product with all 0s, just to start producing s…

imap_processing/cli.py

@@ -64,7 +64,6 @@
 from imap_processing.idex.idex_l1b import idex_l1b
 from imap_processing.idex.idex_l2a import idex_l2a
 from imap_processing.idex.idex_l2b import idex_l2b
-from imap_processing.lo.constants import LoConstants
 from imap_processing.lo.l1a import lo_l1a
 from imap_processing.lo.l1b import lo_l1b
 from imap_processing.lo.l1c import lo_l1c
@@ -1212,53 +1211,6 @@ def do_processing(
 class Lo(ProcessInstrument):
     """Process IMAP-Lo."""
 
-    def pre_processing(self) -> ProcessingInputCollection:
-        """
-        Complete pre-processing.
-
-        Extends the base pre-processing by filtering Lo PSET science inputs to
-        only those whose ``pivot_angle`` is within
-        ``LoConstants.PSET_PIVOT_ANGLE_TOLERANCE`` of
-        ``LoConstants.PSET_PIVOT_ANGLE``. PSET files that fall outside this
-        range are dropped before processing begins.
-
-        Returns
-        -------
-        dependencies : ProcessingInputCollection
-            Object containing dependencies to process.
-        """
-        datasets = super().pre_processing()
-        new_datasets = ProcessingInputCollection()
-
-        for processing_input in datasets.get_processing_inputs():
-            if (
-                processing_input.source == "lo"
-                and processing_input.descriptor == "pset"
-            ):
-                valid_filenames = []
-                for imap_file_path in processing_input.imap_file_paths:
-                    pset = load_cdf(imap_file_path.construct_path())
-                    if "pivot_angle" in pset:
-                        if (
-                            abs(
-                                pset["pivot_angle"].item()
-                                - LoConstants.PSET_PIVOT_ANGLE
-                            )
-                            < LoConstants.PSET_PIVOT_ANGLE_TOLERANCE
-                        ):
-                            valid_filenames.append(str(imap_file_path.filename))
-                        else:
-                            logger.info(
-                                f"Dropping pset {imap_file_path.filename} because "
-                                f"pivot angle is not in range."
-                            )
-                if valid_filenames:
-                    new_datasets.add(type(processing_input)(*valid_filenames))
-            else:
-                new_datasets.add(processing_input)
-
-        return new_datasets
-
     def do_processing(
         self, dependencies: ProcessingInputCollection
     ) -> list[xr.Dataset]:
@@ -1323,18 +1275,21 @@ def do_processing(
             datasets = lo_l1c.lo_l1c(data_dict, anc_dependencies)
 
         elif self.data_level == "l2":
-            data_dict = {}
+            sci_dependencies: dict[str, list[xr.Dataset]] = {}
             science_files = dependencies.get_file_paths(source="lo", descriptor="pset")
+            science_files += dependencies.get_file_paths(source="lo", data_type="l1b")
             anc_dependencies = dependencies.get_file_paths(data_type="ancillary")
 
-            # Load all pset files into datasets
-            if not science_files:
-                logger.info("No valid psets found for L2 processing.")
-                return datasets
+            # Load every input, at every pivot angle, grouped by product.
+            for file in science_files:
+                dataset = load_cdf(file)
+                sci_dependencies.setdefault(dataset.attrs["Logical_source"], []).append(
+                    dataset
+                )
 
-            psets = [load_cdf(file) for file in science_files]
-            data_dict[psets[0].attrs["Logical_source"]] = psets
-            datasets = lo_l2.lo_l2(data_dict, anc_dependencies, self.descriptor)
+            datasets = lo_l2.lo_l2(
+                sci_dependencies, anc_dependencies, self.descriptor, self.start_date
+            )
         return datasets
 
 

imap_processing/lo/constants.py

@@ -8,11 +8,9 @@
 class LoConstants:
     """Constants for Lo which can be used across different levels."""
 
-    # Expected pivot angle [degrees] for pointing sets for generating map products.
-    PSET_PIVOT_ANGLE: float = 90.0
-    # Absolute tolerance [degrees] for accepting a pset's pivot angle
-    # as sufficiently close to the required value.
-    PSET_PIVOT_ANGLE_TOLERANCE: float = 45.0
+    # Absolute tolerance [degrees] for accepting a pset's pivot angle as
+    # sufficiently close to the pivot angle of the map being made.
+    PSET_PIVOT_ANGLE_TOLERANCE: float = 5.0
 
     # Ion species tracked. "H" is mandatory (and should be the first element);
     # any others for which we have histrates may be added here.

imap_processing/lo/l2/lo_l2.py

@@ -18,9 +18,14 @@
     get_pset_directional_mask,
     interpolate_map_flux_to_helio_frame,
 )
-from imap_processing.ena_maps.utils.naming import MapDescriptor
+from imap_processing.ena_maps.utils.naming import DAYS_IN_MONTH, MapDescriptor
 from imap_processing.lo import lo_ancillary
-from imap_processing.spice.time import et_to_datetime64, ttj2000ns_to_et
+from imap_processing.lo.constants import LoConstants
+from imap_processing.spice.time import (
+    et_to_datetime64,
+    str_yyyymmdd_to_ttj2000ns,
+    ttj2000ns_to_et,
+)
 
 logger = logging.getLogger(__name__)
 
@@ -30,48 +35,67 @@
 
 
 def lo_l2(
-    sci_dependencies: dict, anc_dependencies: list, descriptor: str
+    sci_dependencies: dict,
+    anc_dependencies: list,
+    descriptor: str,
+    start_date: str | None = None,
 ) -> list[xr.Dataset]:
     """
     Process IMAP-Lo L1C data into L2 CDF data products.
 
     This is the main entry point for L2 processing. It orchestrates the entire
     processing pipeline from L1C pointing sets to L2 sky maps with intensities.
 
+    Only the pointing sets taken at the pivot angle of the map being made are
+    projected onto it, see ``_inputs_at_map_pivot_angle``. When none of the
+    psets match (or no psets were passed-in to begin with), a quickmap is
+    produced instead.
+
     Parameters
     ----------
     sci_dependencies : dict
         Dictionary of datasets needed for L2 data product creation in xarray Datasets.
-        Must contain "imap_lo_l1c_pset" key with list of pointing set datasets.
+        Should contain "imap_lo_l1c_pset" key with list of pointing set datasets,
+        at any pivot angle.
     anc_dependencies : list
         List of ancillary file paths needed for L2 data product creation.
         Should include efficiency factor files.
     descriptor : str
         The map descriptor to be produced
         (e.g., "ilo90-ena-h-sf-nsp-full-hae-6deg-3mo").
+    start_date : str, optional
+        The start of the map window in YYYYMMDD format. Only required for
+        quickmaps, which have no pointing sets to get their time coverage from.
 
     Returns
     -------
     list[xr.Dataset]
         List containing the processed L2 dataset with rates, intensities,
-        and uncertainties.
+        and uncertainties, or the quickmap if there were no pointing sets at
+        the pivot angle of the map.
 
     Raises
     ------
     ValueError
-        If no pointing set data found in science dependencies.
+        If a quickmap is required but no start_date was given.
     NotImplementedError
         If HEALPix map output is requested (only rectangular maps supported).
     """
     logger.info("Starting IMAP-Lo L2 processing pipeline")
-    if "imap_lo_l1c_pset" not in sci_dependencies:
-        raise ValueError("No pointing set data found in science dependencies")
-    psets = sci_dependencies["imap_lo_l1c_pset"]
 
     # Parse the map descriptor to get species and other attributes
     map_descriptor = MapDescriptor.from_string(descriptor)
     logger.info(f"Processing map for species: {map_descriptor.species}")
 
+    psets = _inputs_at_map_pivot_angle(
+        sci_dependencies.get("imap_lo_l1c_pset", []), map_descriptor
+    )
+    if not psets:
+        logger.info("No psets - trying to create a quickmap.")
+        if start_date is None:
+            raise ValueError(f"A start_date is required to create the map {descriptor}")
+        return _lo_l2_quickmap(sci_dependencies, descriptor, start_date)
+
     # Determine if corrections are needed and prepare oxygen data if required
     (
         sputtering_correction,
@@ -116,6 +140,161 @@ def lo_l2(
     return [dataset]
 
 
+def _inputs_at_map_pivot_angle(
+    datasets: list[xr.Dataset], map_descriptor: MapDescriptor
+) -> list[xr.Dataset]:
+    """
+    Keep only the inputs taken at the pivot angle the map is for.
+
+    Every input we have for the map window is passed in, whatever pivot
+    angle it was taken at. A map is made from one pivot angle only, which is
+    the sensor field of its descriptor ("l090" is the 90 degree pivot angle).
+    Inputs within ``LoConstants.PSET_PIVOT_ANGLE_TOLERANCE`` of that angle are
+    kept, and inputs with no pivot angle at all are dropped.
+
+    Parameters
+    ----------
+    datasets : list[xr.Dataset]
+        The input products available for the map window. Every Lo product
+        records the pivot angle it was taken at.
+    map_descriptor : MapDescriptor
+        The parsed descriptor of the map being made.
+
+    Returns
+    -------
+    list[xr.Dataset]
+        The inputs that belong on this map.
+    """
+    if not isinstance(map_descriptor.sensor, int):
+        # No pivot angle in the descriptor to select inputs with
+        return datasets
+
+    kept = []
+    for dataset in datasets:
+        if "pivot_angle" not in dataset:
+            logger.info("Dropping input with no pivot angle.")
+            continue
+        pivot_angle = dataset["pivot_angle"].item()
+        if (
+            abs(pivot_angle - map_descriptor.sensor)
+            < LoConstants.PSET_PIVOT_ANGLE_TOLERANCE
+        ):
+            kept.append(dataset)
+        else:
+            logger.info(f"Dropping input with pivot angle {pivot_angle}")
+
+    return kept
+
+
+def _lo_l2_quickmap(
+    sci_dependencies: dict, descriptor: str, start_date: str
+) -> list[xr.Dataset]:
+    """
+    Create a correctly shaped, zero-filled L2 quickmap.
+
+    Parameters
+    ----------
+    sci_dependencies : dict
+        Dictionary of the input datasets, keyed by logical source, at any
+        pivot angle.
+    descriptor : str
+        The map descriptor to be produced
+        (e.g., "l090-enansnbs-h-sf-nsp-ram-hae-6deg-6mo").
+    start_date : str
+        The start of the map window in YYYYMMDD format.
+
+    Returns
+    -------
+    list[xr.Dataset]
+        List containing the single quickmap dataset.
+
+    Raises
+    ------
+    NotImplementedError
+        If a HEALPix map is requested (only rectangular maps supported for Lo),
+        or if the map is of a quantity whose variables are not defined yet.
+    """
+    logger.info(f"Creating IMAP-Lo L2 quickmap for descriptor: {descriptor}")
+    map_descriptor = MapDescriptor.from_string(descriptor)
+
+    sky_map = map_descriptor.to_empty_map()
+    if not isinstance(sky_map, RectangularSkyMap):
+        raise NotImplementedError("HEALPix map output not supported for Lo")
+
+    # TODO: No CDF variable attributes are defined for the variables of the
+    # quantities in non-ENA maps, e.g. isn_rate for an ISN map.
+    if not map_descriptor.principal_data.startswith("ena"):
+        raise NotImplementedError(
+            f"Cannot make a quickmap of {map_descriptor.principal_data} for "
+            f"{descriptor}. No CDF variable attributes are defined for "
+            f"{map_descriptor.principal_data_var} in the imap_enamaps_l2 "
+            f"variable attribute files."
+        )
+
+    for logical_source, datasets in sci_dependencies.items():
+        selected = _inputs_at_map_pivot_angle(datasets, map_descriptor)
+        logger.info(
+            f"{len(selected)} of {len(datasets)} {logical_source} inputs are at "
+            f"the pivot angle of {descriptor}"
+        )
+
+    # No pointing sets available, so the epoch bounds are set using `start_date`.
+    duration = cast(str, map_descriptor.duration)
+    if duration.endswith("yr"):
+        duration_months = int(duration.removesuffix("yr")) * 12
+    else:
+        duration_months = int(duration.removesuffix("mo"))
+    sky_map.min_epoch = int(str_yyyymmdd_to_ttj2000ns(start_date))
+    sky_map.max_epoch = sky_map.min_epoch + int(
+        np.timedelta64(duration_months * DAYS_IN_MONTH, "D") / np.timedelta64(1, "ns")
+    )
+
+    # TODO: Figure out how to handle esa_mode properly
+    energies = reduce_geometric_factor_dataset(map_descriptor.species, esa_mode=0)[
+        "Cntr_E"
+    ].values
+
+    # The variables an ENA map has, all with dims (epoch, energy, pixel).
+    float_variables = (
+        "ena_intensity",
+        "ena_intensity_stat_uncert",
+        "ena_intensity_sys_err",
+        "ena_count_rate",
+        "ena_count_rate_stat_uncert",
+        "ena_count",
+        "bg_rate",
+        "bg_rate_stat_uncert",
+        "bg_rate_sys_err",
+        "bg_intensity",
+        "bg_intensity_stat_uncert",
+        "bg_intensity_sys_err",
+        "exposure_factor",
+    )
+    int_variables = ("obs_date", "obs_date_range")
+    variable_dtypes = dict.fromkeys(float_variables, np.float32) | dict.fromkeys(
+        int_variables, np.int64
+    )
+
+    for variable, dtype in variable_dtypes.items():
+        sky_map.data_1d[variable] = xr.DataArray(
+            np.zeros((1, energies.size, sky_map.num_points), dtype=dtype),
+            dims=["epoch", "energy", "pixel"],
+            coords={"energy": energies},
+        )
+
+    dataset = add_geometric_factors(sky_map.to_dataset(), map_descriptor.species)
+    dataset = cleanup_intermediate_variables(dataset)
+
+    return [
+        sky_map.build_cdf_dataset(
+            instrument="lo",
+            level="l2",
+            descriptor=descriptor,
+            external_map_dataset=dataset,
+        )
+    ]
+
+
 def _prepare_corrections(
     map_descriptor: MapDescriptor,
     descriptor: str,