multiBand.py

Go to the documentation of this file.

# This file is part of pipe_tasks.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
__all__ = ["DetectCoaddSourcesConfig", "DetectCoaddSourcesTask",
           "MeasureMergedCoaddSourcesConfig", "MeasureMergedCoaddSourcesTask",
           "DEEP_COADD_BACKGROUND_DOCSTRING",
           ]
 
import dataclasses
 
import astropy.units
import numpy as np
 
from lsst.geom import Extent2I
from lsst.pipe.base import (
    AnnotatedPartialOutputsError,
    Struct,
    PipelineTask,
    PipelineTaskConfig,
    PipelineTaskConnections
)
import lsst.pipe.base.connectionTypes as cT
from lsst.pex.config import Field, ChoiceField, ConfigurableField
from lsst.cell_coadds import MultipleCellCoadd
from lsst.images import Mask, get_legacy_deep_coadd_mask_planes
from lsst.images.cells import CellCoadd
from lsst.images.fields import field_from_legacy_background
from lsst.meas.algorithms import (
    DynamicDetectionTask,
    ExceedsMaxVarianceScaleError,
    InsufficientSourcesError,
    PsfGenerationError,
    ScaleVarianceTask,
    SetPrimaryFlagsTask,
    TooManyMaskedPixelsError,
    ZeroFootprintError,
)
from lsst.meas.base import (
    SingleFrameMeasurementTask,
    ApplyApCorrTask,
    CatalogCalculationTask,
    SkyMapIdGeneratorConfig,
)
from lsst.meas.extensions.scarlet.io import updateCatalogFootprints
from lsst.pipe.tasks.propagateSourceFlags import PropagateSourceFlagsTask
import lsst.afw.image as afwImage
import lsst.afw.math as afwMath
import lsst.afw.table as afwTable
from lsst.daf.base import PropertyList
from lsst.skymap import BaseSkyMap
 
# NOTE: these imports are a convenience so multiband users only have to import this file.
from .mergeDetections import MergeDetectionsConfig, MergeDetectionsTask  # noqa: F401
from .mergeMeasurements import MergeMeasurementsConfig, MergeMeasurementsTask  # noqa: F401
from .multiBandUtils import CullPeaksConfig  # noqa: F401
from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiConfig  # noqa: F401
from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiTask  # noqa: F401
 
 
"""
New set types:
* deepCoadd_det: detections from what used to be processCoadd (tract, patch, filter)
* deepCoadd_mergeDet: merged detections (tract, patch)
* deepCoadd_meas: measurements of merged detections (tract, patch, filter)
* deepCoadd_ref: reference sources (tract, patch)
All of these have associated *_schema catalogs that require no data ID and hold no records.
 
In addition, we have a schema-only dataset, which saves the schema for the PeakRecords in
the mergeDet, meas, and ref dataset Footprints:
* deepCoadd_peak_schema
"""
 
 

 
# The default for DetectCoaddSourcesConfig.backgroundDescription:
DEEP_COADD_BACKGROUND_DOCSTRING = (
    "Background subtracted from the image when generating the Object catalog. "
    "This intentionally oversubtracts the background to reduce blending and ensure "
    "scattered light is subtracted. "
    "Restoring this background does not restore all original backgrounds, "
    "as the coadd was built from background-subtracted visit images; in most "
    "cases this background term is actually quite small "
)
 
 
class DetectCoaddSourcesConnections(PipelineTaskConnections,
                                    dimensions=("tract", "patch", "band", "skymap"),
                                    defaultTemplates={"inputCoaddName": "deep", "outputCoaddName": "deep"}):
    detectionSchema = cT.InitOutput(
        doc="Schema of the detection catalog",
        name="{outputCoaddName}Coadd_det_schema",
        storageClass="SourceCatalog",
    )
    exposure = cT.Input(
        doc="Exposure on which detections are to be performed (if useCellCoadd=False). ",
        name="{inputCoaddName}Coadd",
        storageClass="ExposureF",
        dimensions=("tract", "patch", "band", "skymap")
    )
    exposure_cells = cT.Input(
        doc="Exposure on which detections are to be performed (if useCellCoadd=True). ",
        name="{inputCoaddName}CoaddCell",
        storageClass="MultipleCellCoadd",
        dimensions=("tract", "patch", "band", "skymap"),
    )
    skyMap = cT.Input(
        doc="Description of the skymap's tracts and patches.",
        name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
        storageClass="SkyMap",
        dimensions=("skymap",),
    )
    outputBackgrounds = cT.Output(
        doc="Output Backgrounds used in detection",
        name="{outputCoaddName}Coadd_calexp_background",
        storageClass="Background",
        dimensions=("tract", "patch", "band", "skymap")
    )
    outputSources = cT.Output(
        doc="Detected sources catalog",
        name="{outputCoaddName}Coadd_det",
        storageClass="SourceCatalog",
        dimensions=("tract", "patch", "band", "skymap")
    )
    outputExposure = cT.Output(
        doc="Exposure post detection",
        name="{outputCoaddName}Coadd_calexp",
        storageClass="ExposureF",
        dimensions=("tract", "patch", "band", "skymap")
    )
 
    def __init__(self, *, config=None):
        super().__init__(config=config)
        assert isinstance(config, DetectCoaddSourcesConfig)
 
        if config.imageType == "future":
            if config.useCellCoadds:
                self.exposure_cells = dataclasses.replace(self.exposure_cells, storageClass="CellCoadd")
            else:
                self.exposure = dataclasses.replace(self.exposure, storageClass="CellCoadd")
            self.outputExposure = dataclasses.replace(self.outputExposure, storageClass="CellCoadd")
            del self.outputBackgrounds
        if config.useCellCoadds:
            del self.exposure
        else:
            del self.exposure_cells
        if not self.config.forceExactBinning:
            del self.skyMap
        if self.config.writeOnlyBackgrounds:
            del self.outputExposure
            del self.outputSources
            del self.detectionSchema
 
 
class DetectCoaddSourcesConfig(PipelineTaskConfig, pipelineConnections=DetectCoaddSourcesConnections):
    """Configuration parameters for the DetectCoaddSourcesTask
    """
 
    doScaleVariance = Field(dtype=bool, default=True, doc="Scale variance plane using empirical noise?")
    scaleVariance = ConfigurableField(target=ScaleVarianceTask, doc="Variance rescaling")
    detection = ConfigurableField(target=DynamicDetectionTask, doc="Source detection")
    coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
    useCellCoadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?")
    hasFakes = Field(
        dtype=bool,
        default=False,
        doc="Should be set to True if fake sources have been inserted into the input data.",
    )
    idGenerator = SkyMapIdGeneratorConfig.make_field()
    forceExactBinning = Field(
        dtype=bool,
        default=False,
        doc=(
            "Check that the background bin size evenly divides the patch inner region, and "
            "crop the outer region to an integer number of bins."
        )
    )
    writeOnlyBackgrounds = Field(dtype=bool, default=False, doc="If true, only save the background models.")
    writeEmptyBackgrounds = Field(
        dtype=bool,
        default=True,
        doc=(
            "If true, save a placeholder background with NaNs in all bins (but the right geometry) when "
            "there are no pixels to compute a background from.  This can be useful if a later task combines "
            "backgrounds from multiple patches as input."
        )
    )
    imageType = ChoiceField(
        "Which image type to use for the input and output coadd. "
        "This option only directly affects connection storage classes and hence 'runQuantum'; the 'run' "
        "method behavior is determined by which type is actually passed in.",
        allowed={
            "legacy": (
                "Read a lsst.cell_coadds.MultipleCellCoadd (if useCellCoadd) or "
                "lsst.afw.image.Exposure (if not useCellCoadd) and write an lsst.afw.image.Exposure."
            ),
            "future": (
                "Read and write lsst.images.cells.CellCoadd (useCellCoadd just "
                "sets which of 'exposure' or 'exposure_cells' will be used for inputs). "
                "The 'outputBackground' is deleted, writeEmptyBackgrounds is ignored, and "
                "writeOnlyBackgrounds=True is invalid."
            ),
        },
        dtype=str,
        optional=False,
        default="legacy",
    )
    backgroundName = Field(
        "Name of the subtracted background, to be stored with the image when the input and "
        "output images are lsst.images.cells.CellCoadd.",
        dtype=str,
        default="object",
    )
    backgroundDescription = Field(
        "Description of the subtracted background, to be stored with the image when the input and "
        "output images are lsst.images.cells.CellCoadd.",
        dtype=str,
        default=DEEP_COADD_BACKGROUND_DOCSTRING,
    )
 
    def setDefaults(self):
        super().setDefaults()
        self.detection.thresholdType = "pixel_stdev"
        self.detection.isotropicGrow = True
        # Coadds are made from background-subtracted CCDs, so any background subtraction should be very basic
        self.detection.reEstimateBackground = False
        self.detection.background.useApprox = False
        self.detection.background.binSize = 4096
        self.detection.background.undersampleStyle = 'REDUCE_INTERP_ORDER'
        self.detection.doTempWideBackground = True  # Suppress large footprints that overwhelm the deblender
        # Include band in packed data IDs that go into object IDs (None -> "as
        # many bands as are defined", rather than the default of zero).
        self.idGenerator.packer.n_bands = None
 
    def validate(self):
        super().validate()
        if self.imageType == "future":
            if self.doScaleVariance:
                raise ValueError("doScaleVariance=True is not compatible with imageType='future'")
            if self.forceExactBinning:
                raise ValueError("forceExactBinning=True is not compatible with imageType='future'")
            if self.writeOnlyBackgrounds:
                raise ValueError("writeOnlyBackgrounds=True is not compatible with imageType='future'")
 
 
class DetectCoaddSourcesTask(PipelineTask):
    """Detect sources on a single filter coadd.
 
    Coadding individual visits requires each exposure to be warped. This
    introduces covariance in the noise properties across pixels. Before
    detection, we correct the coadd variance by scaling the variance plane in
    the coadd to match the observed variance. This is an approximate
    approach -- strictly, we should propagate the full covariance matrix --
    but it is simple and works well in practice.
 
    After scaling the variance plane, we detect sources and generate footprints
    by delegating to the @ref SourceDetectionTask_ "detection" subtask.
 
    DetectCoaddSourcesTask is meant to be run after assembling a coadded image
    in a given band. The purpose of the task is to update the background,
    detect all sources in a single band and generate a set of parent
    footprints. Subsequent tasks in the multi-band processing procedure will
    merge sources across bands and, eventually, perform forced photometry.
 
    Parameters
    ----------
    schema : `lsst.afw.table.Schema`, optional
        Initial schema for the output catalog, modified-in place to include all
        fields set by this task.  If None, the source minimal schema will be used.
    **kwargs
        Additional keyword arguments.
    """
 
    _DefaultName = "detectCoaddSources"
    ConfigClass = DetectCoaddSourcesConfig
 
    def __init__(self, schema=None, **kwargs):
        # N.B. Super is used here to handle the multiple inheritance of PipelineTasks, the init tree
        # call structure has been reviewed carefully to be sure super will work as intended.
        super().__init__(**kwargs)
        if schema is None:
            schema = afwTable.SourceTable.makeMinimalSchema()
        self.schema = schema
        self.makeSubtask("detection", schema=self.schema)
        if self.config.doScaleVariance:
            self.makeSubtask("scaleVariance")
 
        self.detectionSchema = afwTable.SourceCatalog(self.schema)
 
    def runQuantum(self, butlerQC, inputRefs, outputRefs):
        inputs = butlerQC.get(inputRefs)
        idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
 
        if self.config.useCellCoadds:
            multiple_cell_coadd = inputs.pop("exposure_cells")
            match self.config.imageType:
                case "legacy":
                    exposure = multiple_cell_coadd.stitch().asExposure()
                case "future":
                    exposure = multiple_cell_coadd  # conversion deferred to run().
                case _:
                    raise AssertionError(f"Invalid choice {self.config.imageType!r} for imageType.")
        else:
            exposure = inputs.pop("exposure")
 
        skyMap = inputs.pop("skyMap", None)
        if skyMap is not None:
            patchInfo = skyMap[butlerQC.quantum.dataId["tract"]][butlerQC.quantum.dataId["patch"]]
        else:
            patchInfo = None
 
        assert not inputs, "runQuantum got more inputs than expected."
        try:
            outputs = self.run(
                exposure=exposure,
                idFactory=idGenerator.make_table_id_factory(),
                expId=idGenerator.catalog_id,
                patchInfo=patchInfo,
            )
        except (
            TooManyMaskedPixelsError,
            ExceedsMaxVarianceScaleError,
            InsufficientSourcesError,
            PsfGenerationError,
            ZeroFootprintError,
        ) as e:
            if self.config.writeEmptyBackgrounds:
                butlerQC.put(self._makeEmptyBackground(exposure, patchInfo), outputRefs.outputBackgrounds)
            # Detection failed, so clear any leftover detected mask planes.
            if isinstance(exposure, CellCoadd):
                # If we passed a CellCoadd in, it won't be modified until 'run'
                # is about to exit, and hence it can't have picked up a
                # DETECTED_NEGATIVE plane, which can only be temporary here.
                # But it might have a preexisting DETECTED plane (i.e. a union
                # of the DETECTED plane from the warps) that we'd still want to
                # clear.
                exposure.mask.clear("DETECTED")
            else:
                for maskName in ["DETECTED", "DETECTED_NEGATIVE"]:
                    if maskName in exposure.mask.getMaskPlaneDict().keys():
                        detectedMask = exposure.mask.getMaskPlane(maskName)
                        exposure.mask.clearMaskPlane(detectedMask)
            error = AnnotatedPartialOutputsError.annotate(
                e,
                self,
                exposure,
                log=self.log,
            )
            butlerQC.put(exposure, outputRefs.outputExposure)
            raise error from e
 
        butlerQC.put(outputs, outputRefs)
 
    def run(self, exposure, idFactory, expId, patchInfo=None):
        """Run detection on an exposure.
 
        First scale the variance plane to match the observed variance
        using ``ScaleVarianceTask``. Then invoke the ``SourceDetectionTask_`` "detection" subtask to
        detect sources.
 
        Parameters
        ----------
        exposure : `lsst.afw.image.Exposure` or `lsst.images.cells.CellCoadd`.
            Exposure on which to detect (may be background-subtracted and scaled,
            depending on configuration).
        idFactory : `lsst.afw.table.IdFactory`
            IdFactory to set source identifiers.
        expId : `int`
            Exposure identifier (integer) for RNG seed.
        patchInfo : `lsst.skymap.PatchInfo`, optional
            Description of the patch geometry.  Only needed if
            `~DetectCoaddSourceConfig.forceExactBinning` is `True`.
 
        Returns
        -------
        result : `lsst.pipe.base.Struct`
            Results as a struct with attributes:
 
            ``outputSources``
                Catalog of detections (`lsst.afw.table.SourceCatalog`).
            ``outputBackgrounds``
                List of backgrounds (`list`).
            ``outputExposure``
                The background-subtracted coadd image, with its mask plane
                updated to include detections.  This will have the same type
                as ``exposure``.
        """
        cell_coadd = None
        if isinstance(exposure, CellCoadd):
            cell_coadd = exposure
            exposure = cell_coadd.to_legacy()
        if self.config.forceExactBinning:
            if cell_coadd is not None:
                raise ValueError("forceExactBinning=True is not compatible with CellCoadd inputs")
            exposure = self._cropToExactBinning(exposure, patchInfo)
        if self.config.doScaleVariance:
            if cell_coadd is not None:
                raise ValueError("doScaleVariance=True is not compatible with CellCoadd inputs")
            varScale = self.scaleVariance.run(exposure.maskedImage)
            exposure.getMetadata().add("VARIANCE_SCALE", varScale)
        backgrounds = afwMath.BackgroundList()
        table = afwTable.SourceTable.make(self.schema, idFactory)
        detections = self.detection.run(table, exposure, expId=expId)
        sources = detections.sources
        if hasattr(detections, "background") and detections.background:
            for bg in detections.background:
                backgrounds.append(bg)
        if len(backgrounds) == 0:
            # Persist a constant background with value of NaN to get around
            # inability to persist empty BackgroundList.
            emptyBg = self._makeEmptyBackground(exposure, patchInfo)
            backgrounds.append(emptyBg)
        if cell_coadd is not None:
            cell_coadd.image.array[...] = exposure.image.array
            cell_coadd.mask = Mask.from_legacy(
                exposure.mask, plane_map=get_legacy_deep_coadd_mask_planes()
            ).view(sky_projection=cell_coadd.sky_projection)
            if backgrounds:
                cell_coadd.backgrounds._backgrounds.clear()
                cell_coadd.backgrounds.add(
                    self.config.backgroundName,
                    field_from_legacy_background(backgrounds, unit=astropy.units.nJy),
                    self.config.backgroundDescription,
                    is_subtracted=True,
                )
            exposure = cell_coadd
        return Struct(outputSources=sources, outputBackgrounds=backgrounds, outputExposure=exposure)
 
    def _cropToExactBinning(self, exposure, patchInfo):
        """Crop a coadd `~lsst.afw.image.Exposure` instance to ensure exact
        background binning.
 
        Parameters
        ----------
        exposure : `lsst.afw.image.Exposure`
            Exposure to crop, assumed to cover the patch outer bounding box.
        patchInfo : `lsst.skymap.PatchInfo`
            Description of the patch geometry.
 
        Returns
        -------
        cropped : `lsst.afw.image.Exposure`
            View of ``exposure`` with background bins that evenly divide both
            the full cropped image and the patch inner region.  The bounding
            box is guaranteed to contain the patch inner bounding box and be
            contained by the patch outer bounding box.
 
        Raises
        ------
        ValueError
            Raised if the patch inner region width or height is not a multiple
            of the background bin size.
        """
        bbox = patchInfo.getInnerBBox()
        if bbox.width % self.detection.background.binSizeX:
            raise ValueError(
                f"Patch inner width {bbox.width} does not evenly "
                f"divide bin width {self.detection.background.binSizeX}."
            )
        if bbox.height % self.detection.background.binSizeY:
            raise ValueError(
                f"Patch inner height {bbox.height} does not evenly "
                f"divide bin height {self.detection.background.binSizeY}."
            )
        outer_bbox = patchInfo.getOuterBBox()
        n_bins_grow_x = (bbox.x.begin - outer_bbox.x.begin) // self.detection.background.binSizeX
        n_bins_grow_y = (bbox.y.begin - outer_bbox.y.begin) // self.detection.background.binSizeY
        bbox.grow(
            Extent2I(
                n_bins_grow_x*self.detection.background.binSizeX,
                n_bins_grow_y*self.detection.background.binSizeY,
            )
        )
        assert outer_bbox.contains(bbox)
        assert bbox.contains(patchInfo.getInnerBBox())
        assert bbox.width % self.detection.background.binSizeX == 0
        assert bbox.height % self.detection.background.binSizeY == 0
        return exposure[bbox]
 
    def _makeEmptyBackground(self, exposure, patchInfo=None):
        """Construct an empty `lsst.afw.math.BackgroundList` with NaN values.
 
        Parameters
        ----------
        exposure : `lsst.afw.image.Exposure`
            Exposure that the background should correspond to.
        patchInfo : `lsst.skymap.PatchInfo`, optional
            Description of the patch geometry.  Only needed if
            `~DetectCoaddSourceConfig.forceExactBinning` is `True`.
 
        Returns
        -------
        background : `lsst.afw.math.BackgroundList`
            A background object with a single layer and the same bin geometry
            that a background for that exposure would have had if it had enough
            usable pixels.  This object cannot actually be used for background
            subtraction.
        """
        # Create a backgroundList with one entry whose "stats image" is NaNs
        # and has all pixels set as NO_DATA.
        if self.config.forceExactBinning:
            exposure = self._cropToExactBinning(exposure, patchInfo).clone()
 
        bgLevel = np.nan
        bgStats = afwImage.MaskedImageF(1, 1)
        bgStats.set(bgLevel, 0, bgLevel)
        bg = afwMath.BackgroundMI(exposure.getBBox(), bgStats)
        bgData = (bg, afwMath.Interpolate.LINEAR, afwMath.REDUCE_INTERP_ORDER,
                  afwMath.ApproximateControl.UNKNOWN, 0, 0, False)
        background = afwMath.BackgroundList()
        background.append(bgData)
        for bg, *_ in background:
            stats = bg.getStatsImage()
            stats.mask.array[:, :] = stats.mask.getPlaneBitMask("NO_DATA")
            stats.variance.array[:, :] = 0.0
        return background
 
 
class MeasureMergedCoaddSourcesConnections(
    PipelineTaskConnections,
    dimensions=("tract", "patch", "band", "skymap"),
    defaultTemplates={
        "inputCoaddName": "deep",
        "outputCoaddName": "deep",
    },
):
    inputSchema = cT.InitInput(
        doc="Input schema for measure merged task produced by a deblender or detection task",
        name="{inputCoaddName}Coadd_deblendedFlux_schema",
        storageClass="SourceCatalog"
    )
    outputSchema = cT.InitOutput(
        doc="Output schema after all new fields are added by task",
        name="{inputCoaddName}Coadd_meas_schema",
        storageClass="SourceCatalog"
    )
    exposure = cT.Input(
        doc="Input non-cell-based coadd image",
        name="{inputCoaddName}Coadd_calexp",
        storageClass="ExposureF",
        dimensions=("tract", "patch", "band", "skymap")
    )
    exposure_cells = cT.Input(
        doc="Input cell-based coadd image",
        name="{inputCoaddName}CoaddCell",
        storageClass="MultipleCellCoadd",
        dimensions=("tract", "patch", "band", "skymap"),
    )
    background = cT.Input(
        doc="Background to subtract from cell-based coadd image",
        name="{inputCoaddName}Coadd_calexp_background",
        storageClass="Background",
        dimensions=("tract", "patch", "band", "skymap")
    )
    skyMap = cT.Input(
        doc="SkyMap to use in processing",
        name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
        storageClass="SkyMap",
        dimensions=("skymap",),
    )
    sourceTableHandles = cT.Input(
        doc=("Source tables that are derived from the ``CalibrateTask`` sources. "
             "These tables contain astrometry and photometry flags, and optionally "
             "PSF flags."),
        name="sourceTable_visit",
        storageClass="ArrowAstropy",
        dimensions=("instrument", "visit"),
        multiple=True,
        deferLoad=True,
    )
    finalizedSourceTableHandles = cT.Input(
        doc=("Finalized source tables from ``FinalizeCalibrationTask``. These "
             "tables contain PSF flags from the finalized PSF estimation."),
        name="finalized_src_table",
        storageClass="ArrowAstropy",
        dimensions=("instrument", "visit"),
        multiple=True,
        deferLoad=True,
    )
    finalVisitSummaryHandles = cT.Input(
        doc="Final visit summary table",
        name="finalVisitSummary",
        storageClass="ExposureCatalog",
        dimensions=("instrument", "visit"),
        multiple=True,
        deferLoad=True,
    )
    scarletCatalog = cT.Input(
        doc="Catalogs produced by multiband deblending",
        name="{inputCoaddName}Coadd_deblendedCatalog",
        storageClass="SourceCatalog",
        dimensions=("tract", "patch", "skymap"),
    )
    scarletModels = cT.Input(
        doc="Multiband scarlet models produced by the deblender",
        name="{inputCoaddName}Coadd_scarletModelData",
        storageClass="LsstScarletModelData",
        dimensions=("tract", "patch", "skymap"),
    )
    outputSources = cT.Output(
        doc="Source catalog containing all the measurement information generated in this task",
        name="{outputCoaddName}Coadd_meas",
        dimensions=("tract", "patch", "band", "skymap"),
        storageClass="SourceCatalog",
    )
 
    def __init__(self, *, config=None):
        super().__init__(config=config)
        if not config.doPropagateFlags:
            del self.sourceTableHandles
            del self.finalizedSourceTableHandles
            del self.finalVisitSummaryHandles
        else:
            # Check for types of flags required.
            if not config.propagateFlags.source_flags:
                del self.sourceTableHandles
            if not config.propagateFlags.finalized_source_flags:
                del self.finalizedSourceTableHandles
        if not config.doAddFootprints:
            del self.scarletModels
        if self.config.imageType == "future":
            self.exposure = dataclasses.replace(self.exposure, storageClass="CellCoadd")
            del self.exposure_cells
            del self.background
        elif self.config.useCellCoadds:
            del self.exposure
        else:
            del self.exposure_cells
            del self.background
 
 
class MeasureMergedCoaddSourcesConfig(PipelineTaskConfig,
                                      pipelineConnections=MeasureMergedCoaddSourcesConnections):
    """Configuration parameters for the MeasureMergedCoaddSourcesTask
    """
    doAddFootprints = Field(dtype=bool,
                            default=True,
                            doc="Whether or not to add footprints to the input catalog from scarlet models. "
                                "This should be true whenever using the multi-band deblender, "
                                "otherwise this should be False.")
    doConserveFlux = Field(dtype=bool, default=True,
                           doc="Whether to use the deblender models as templates to re-distribute the flux "
                               "from the 'exposure' (True), or to perform measurements on the deblender "
                               "model footprints.")
    doStripFootprints = Field(dtype=bool, default=True,
                              doc="Whether to strip footprints from the output catalog before "
                                  "saving to disk. "
                                  "This is usually done when using scarlet models to save disk space.")
    useCellCoadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?")
    measurement = ConfigurableField(target=SingleFrameMeasurementTask, doc="Source measurement")
    setPrimaryFlags = ConfigurableField(target=SetPrimaryFlagsTask, doc="Set flags for primary tract/patch")
    doPropagateFlags = Field(
        dtype=bool, default=True,
        doc="Whether to match sources to CCD catalogs to propagate flags (to e.g. identify PSF stars)"
    )
    propagateFlags = ConfigurableField(target=PropagateSourceFlagsTask, doc="Propagate source flags to coadd")
    coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
    psfCache = Field(dtype=int, default=100, doc="Size of psfCache")
    checkUnitsParseStrict = Field(
        doc="Strictness of Astropy unit compatibility check, can be 'raise', 'warn' or 'silent'",
        dtype=str,
        default="raise",
    )
    doApCorr = Field(
        dtype=bool,
        default=True,
        doc="Apply aperture corrections"
    )
    applyApCorr = ConfigurableField(
        target=ApplyApCorrTask,
        doc="Subtask to apply aperture corrections"
    )
    doRunCatalogCalculation = Field(
        dtype=bool,
        default=True,
        doc='Run catalogCalculation task'
    )
    catalogCalculation = ConfigurableField(
        target=CatalogCalculationTask,
        doc="Subtask to run catalogCalculation plugins on catalog"
    )
 
    hasFakes = Field(
        dtype=bool,
        default=False,
        doc="Should be set to True if fake sources have been inserted into the input data."
    )
    idGenerator = SkyMapIdGeneratorConfig.make_field()
    imageType = ChoiceField(
        "Which image type to expect for the input coadd. "
        "This option only directly affects connection storage classes and hence 'runQuantum'; the 'run' "
        "method behavior is determined by which type is actually passed in.",
        allowed={
            "legacy": (
                "Read a lsst.cell_coadds.MultipleCellCoadd via 'exposure_cells` and restore 'background' "
                "(if useCellCoadd) or lsst.afw.image.Exposure via `exposure` (if not useCellCoadd)."
            ),
            "future": (
                "Read lsst.images.cells.CellCoadd via the 'exposure' connection.  useCellCoadd is ignored."
            ),
        },
        dtype=str,
        optional=False,
        default="legacy",
    )
 
    def setDefaults(self):
        super().setDefaults()
        self.measurement.plugins.names |= ['base_InputCount',
                                           'base_Variance',
                                           'base_LocalPhotoCalib',
                                           'base_LocalWcs']
 
        # TODO: Remove STREAK in DM-44658, streak masking to happen only in
        # ip_diffim; if we can propagate the streak mask from diffim, we can
        # still set flags with it here.
        self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['CLIPPED', 'SENSOR_EDGE',
                                                                       'INEXACT_PSF']
        self.measurement.plugins['base_PixelFlags'].masksFpCenter = ['CLIPPED', 'SENSOR_EDGE',
                                                                     'INEXACT_PSF']
 
 
class MeasureMergedCoaddSourcesTask(PipelineTask):
    """Deblend sources from main catalog in each coadd seperately and measure.
 
    Use peaks and footprints from a master catalog to perform deblending and
    measurement in each coadd.
 
    Given a master input catalog of sources (peaks and footprints) or deblender
    outputs(including a HeavyFootprint in each band), measure each source on
    the coadd. Repeating this procedure with the same master catalog across
    multiple coadds will generate a consistent set of child sources.
 
    The deblender retains all peaks and deblends any missing peaks (dropouts in
    that band) as PSFs. Source properties are measured and the @c is-primary
    flag (indicating sources with no children) is set. Visit flags are
    propagated to the coadd sources.
 
    After MeasureMergedCoaddSourcesTask has been run on multiple coadds, we
    have a set of per-band catalogs. The next stage in the multi-band
    processing procedure will merge these measurements into a suitable catalog
    for driving forced photometry.
 
    Parameters
    ----------
    schema : ``lsst.afw.table.Schema`, optional
        The schema of the merged detection catalog used as input to this one.
    peakSchema : ``lsst.afw.table.Schema`, optional
        The schema of the PeakRecords in the Footprints in the merged detection catalog.
    initInputs : `dict`, optional
        Dictionary that can contain a key ``inputSchema`` containing the
        input schema. If present will override the value of ``schema``.
    **kwargs
        Additional keyword arguments.
    """
 
    _DefaultName = "measureCoaddSources"
    ConfigClass = MeasureMergedCoaddSourcesConfig
 
    def __init__(self, schema=None, peakSchema=None, initInputs=None, **kwargs):
        super().__init__(**kwargs)
        if initInputs is not None:
            schema = initInputs['inputSchema'].schema
        if schema is None:
            raise ValueError("Schema must be defined.")
        self.schemaMapper = afwTable.SchemaMapper(schema)
        self.schemaMapper.addMinimalSchema(schema)
        self.schema = self.schemaMapper.getOutputSchema()
        self.algMetadata = PropertyList()
        self.makeSubtask("measurement", schema=self.schema, algMetadata=self.algMetadata)
        self.makeSubtask("setPrimaryFlags", schema=self.schema)
        if self.config.doPropagateFlags:
            self.makeSubtask("propagateFlags", schema=self.schema)
        self.schema.checkUnits(parse_strict=self.config.checkUnitsParseStrict)
        if self.config.doApCorr:
            self.makeSubtask("applyApCorr", schema=self.schema)
        if self.config.doRunCatalogCalculation:
            self.makeSubtask("catalogCalculation", schema=self.schema)
 
        self.outputSchema = afwTable.SourceCatalog(self.schema)
 
    def runQuantum(self, butlerQC, inputRefs, outputRefs):
        inputs = butlerQC.get(inputRefs)
 
        if self.config.imageType == "future":
            coadd = inputs.pop("exposure")
            band = inputRefs.exposure.dataId["band"]
            # Instead of going directly from lsst.images.cells.CellCoadd to
            # Exposure, it's cleaner for now to go through MultipleCellCoadd
            # because the apCorrMap and ccdInputs need special handling - the
            # cell-based versions can't be attached to Exposure.  Eventually
            # we'll rewrite the lower-level code to use the lsst.images
            # equivalents natively.
            coadd = coadd.to_legacy_cell_coadd()
        elif self.config.useCellCoadds:
            coadd = inputs.pop("exposure_cells")
            band = inputRefs.exposure_cells.dataId["band"]
        else:
            coadd = inputs.pop("exposure")
            band = inputRefs.exposure.dataId["band"]
        if isinstance(coadd, MultipleCellCoadd):
            stitched_coadd = coadd.stitch()
            exposure = stitched_coadd.asExposure()
            if self.config.imageType == "legacy":
                background = inputs.pop("background")
                exposure.image -= background.getImage()
            ccdInputs = stitched_coadd.ccds
            apCorrMap = stitched_coadd.ap_corr_map
        else:
            exposure = coadd
            # Set psfcache only when we don't have a cell-based coadd.
            exposure.getPsf().setCacheCapacity(self.config.psfCache)
 
            ccdInputs = exposure.getInfo().getCoaddInputs().ccds
            apCorrMap = exposure.getInfo().getApCorrMap()
 
        # Get unique integer ID for IdFactory and RNG seeds; only the latter
        # should really be used as the IDs all come from the input catalog.
        idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
 
        # Transform inputCatalog
        table = afwTable.SourceTable.make(self.schema, idGenerator.make_table_id_factory())
        sources = afwTable.SourceCatalog(table)
        # Load the correct input catalog
        inputCatalog = inputs.pop("scarletCatalog")
        catalogRef = inputRefs.scarletCatalog
        sources.extend(inputCatalog, self.schemaMapper)
        del inputCatalog
        # Add the HeavyFootprints to the deblended sources
        if self.config.doAddFootprints:
            modelData = inputs.pop('scarletModels')
            if self.config.doConserveFlux:
                imageForRedistribution = exposure
            else:
                imageForRedistribution = None
            updateCatalogFootprints(
                modelData=modelData,
                catalog=sources,
                band=band,
                imageForRedistribution=imageForRedistribution,
                removeScarletData=True,
                updateFluxColumns=True,
            )
        table = sources.getTable()
        table.setMetadata(self.algMetadata)  # Capture algorithm metadata to write out to the source catalog.
 
        skyMap = inputs.pop('skyMap')
        tractNumber = catalogRef.dataId['tract']
        tractInfo = skyMap[tractNumber]
        patchInfo = tractInfo.getPatchInfo(catalogRef.dataId['patch'])
        skyInfo = Struct(
            skyMap=skyMap,
            tractInfo=tractInfo,
            patchInfo=patchInfo,
            wcs=tractInfo.getWcs(),
            bbox=patchInfo.getOuterBBox()
        )
 
        sourceTableHandleDict = None
        finalizedSourceTableHandleDict = None
        finalVisitSummaryHandleDict = None
        if self.config.doPropagateFlags:
            if "sourceTableHandles" in inputs:
                sourceTableHandles = inputs.pop("sourceTableHandles")
                sourceTableHandleDict = {handle.dataId["visit"]: handle for handle in sourceTableHandles}
            if "finalizedSourceTableHandles" in inputs:
                finalizedSourceTableHandles = inputs.pop("finalizedSourceTableHandles")
                finalizedSourceTableHandleDict = {handle.dataId["visit"]: handle
                                                  for handle in finalizedSourceTableHandles}
            if "finalVisitSummaryHandles" in inputs:
                finalVisitSummaryHandles = inputs.pop("finalVisitSummaryHandles")
                finalVisitSummaryHandleDict = {handle.dataId["visit"]: handle
                                               for handle in finalVisitSummaryHandles}
 
        assert not inputs, "runQuantum got more inputs than expected."
        outputs = self.run(
            exposure=exposure,
            sources=sources,
            skyInfo=skyInfo,
            exposureId=idGenerator.catalog_id,
            ccdInputs=ccdInputs,
            sourceTableHandleDict=sourceTableHandleDict,
            finalizedSourceTableHandleDict=finalizedSourceTableHandleDict,
            finalVisitSummaryHandleDict=finalVisitSummaryHandleDict,
            apCorrMap=apCorrMap,
        )
        # Strip HeavyFootprints to save space on disk
        if self.config.doStripFootprints:
            sources = outputs.outputSources
            for source in sources[sources["parent"] != 0]:
                source.setFootprint(None)
        butlerQC.put(outputs, outputRefs)
 
    def run(self, exposure, sources, skyInfo, exposureId, ccdInputs=None,
            sourceTableHandleDict=None, finalizedSourceTableHandleDict=None, finalVisitSummaryHandleDict=None,
            apCorrMap=None):
        """Run measurement algorithms on the input exposure, and optionally populate the
        resulting catalog with extra information.
 
        Parameters
        ----------
        exposure : `lsst.afw.image.Exposure`
            The input exposure on which measurements are to be performed.
        sources :  `lsst.afw.table.SourceCatalog`
            A catalog built from the results of merged detections, or
            deblender outputs.
        parentCatalog : `lsst.afw.table.SourceCatalog`
            Catalog of parent sources corresponding to sources.
        skyInfo : `lsst.pipe.base.Struct`
            A struct containing information about the position of the input exposure within
            a `SkyMap`, the `SkyMap`, its `Wcs`, and its bounding box.
        exposureId : `int` or `bytes`
            Packed unique number or bytes unique to the input exposure.
        ccdInputs : `lsst.afw.table.ExposureCatalog`, optional
            Catalog containing information on the individual visits which went into making
            the coadd.
        sourceTableHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
            Dict for sourceTable_visit handles (key is visit) for propagating flags.
            These tables contain astrometry and photometry flags, and optionally PSF flags.
        finalizedSourceTableHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
            Dict for finalized_src_table handles (key is visit) for propagating flags.
            These tables contain PSF flags from the finalized PSF estimation.
        finalVisitSummaryHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
            Dict for visit_summary handles (key is visit) for visit-level information.
            These tables contain the WCS information of the single-visit input images.
        apCorrMap : `lsst.afw.image.ApCorrMap`, optional
            Aperture correction map attached to the ``exposure``. If None, it
            will be read from the ``exposure``.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
            Results of running measurement task. Will contain the catalog in the
            sources attribute.
        """
        if self.config.doPropagateFlags:
            # These mask planes may not be defined on the coadds always.
            # We add the mask planes, which is a no-op if already defined.
            for maskPlane in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
                exposure.mask.addMaskPlane(maskPlane)
            for maskPlane in self.config.measurement.plugins["base_PixelFlags"].masksFpCenter:
                exposure.mask.addMaskPlane(maskPlane)
 
        self._ensureMaskPlanes()
 
        self.measurement.run(sources, exposure, exposureId=exposureId)
 
        if self.config.doApCorr:
            if apCorrMap is None:
                apCorrMap = exposure.getInfo().getApCorrMap()
            self.applyApCorr.run(
                catalog=sources,
                apCorrMap=apCorrMap,
            )
 
        # TODO DM-11568: this contiguous check-and-copy could go away if we
        # reserve enough space during SourceDetection and/or SourceDeblend.
        # NOTE: sourceSelectors require contiguous catalogs, so ensure
        # contiguity now, so views are preserved from here on.
        if not sources.isContiguous():
            sources = sources.copy(deep=True)
 
        if self.config.doRunCatalogCalculation:
            self.catalogCalculation.run(sources)
 
        self.setPrimaryFlags.run(sources, skyMap=skyInfo.skyMap, tractInfo=skyInfo.tractInfo,
                                 patchInfo=skyInfo.patchInfo)
        if self.config.doPropagateFlags:
            self.propagateFlags.run(
                sources,
                ccdInputs,
                sourceTableHandleDict,
                finalizedSourceTableHandleDict,
                finalVisitSummaryHandleDict,
            )
 
        results = Struct()
        results.outputSources = sources
        return results
 
    def _ensureMaskPlanes(self):
        """Ensure the global mask dictionary has all of the mask planes
        needed for PixelFlags algorithms.
 
        When mask planes are added, this essentially guarantees that the
        corresponding PixelFlags columns will be wholly False, and usually
        we'd prefer to remove them from the configuration.  But those config
        changes imply a schema changes, and that's not always viable (e.g. on
        a release branch).
        """
        needed = set(self.measurement.plugins["base_PixelFlags"].config.masksFpCenter)
        needed.update(self.measurement.plugins["base_PixelFlags"].config.masksFpAnywhere)
        existing = afwImage.MaskX().getMaskPlaneDict().keys()
        for plane in sorted(needed - existing):
            self.log.warning(
                "Adding mask plane %r with no pixel set to satisfy PixelFlags configuration.", plane
            )
            afwImage.MaskX.addMaskPlane(plane)