Merge pull request #494 from lsst/tickets/DM-29279

DM-29279: Add healsparse input map generation to AssembleCoaddTask

doc/lsst.pipe.tasks/tasks/lsst.pipe.tasks.healSparseMapping.HealSparseInputMapTask.rst

@@ -0,0 +1,26 @@
+.. lsst-task-topic:: lsst.pipe.tasks.healSparseMapping.HealSparseInputMapTask
+
+######################
+HealSparseInputMapTask
+######################
+
+.. _lsst.pipe.tasks.healSparseMapping.HealSparseInputMapTask-api:
+
+Python API summary
+==================
+
+.. lsst-task-api-summary:: lsst.pipe.tasks.healSparseMapping.HealSparseInputMapTask
+
+.. _lsst.pipe.tasks.healSparseMapping.HealSparseInputMapTask-subtasks:
+
+Retargetable subtasks
+=====================
+
+.. lsst-task-config-subtasks:: lsst.pipe.tasks.healSparseMapping.HealSparseInputMapTask
+
+.. _lsst.pipe.tasks.healSparseMapping.HealSparseInputMapTask-fields:
+
+Configuration fields
+====================
+
+.. lsst-task-config-fields:: lsst.pipe.tasks.healSparseMapping.HealSparseInputMapTask

pipelines/DRP.yaml

@@ -16,7 +16,10 @@ tasks:
     config:
       makePsfMatched: true
       doWriteEmptyWarps: true
-  assembleCoadd: lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask
+  assembleCoadd:
+    class: lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask
+    config:
+      doInputMap: true
   detection: lsst.pipe.tasks.multiBand.DetectCoaddSourcesTask
   mergeDetections: lsst.pipe.tasks.mergeDetections.MergeDetectionsTask
   deblend: lsst.pipe.tasks.deblendCoaddSourcesPipeline.DeblendCoaddSourcesMultiTask

python/lsst/pipe/tasks/assembleCoadd.py

@@ -44,6 +44,7 @@
 from .coaddHelpers import groupPatchExposures, getGroupDataRef
 from .scaleVariance import ScaleVarianceTask
 from .maskStreaks import MaskStreaksTask
+from .healSparseMapping import HealSparseInputMapTask
 from lsst.meas.algorithms import SourceDetectionTask
 from lsst.daf.butler import DeferredDatasetHandle
 
@@ -99,6 +100,12 @@ class AssembleCoaddConnections(pipeBase.PipelineTaskConnections,
         storageClass="ImageU",
         dimensions=("tract", "patch", "skymap", "band"),
     )
+    inputMap = pipeBase.connectionTypes.Output(
+        doc="Output healsparse map of input images",
+        name="{outputCoaddName}Coadd_inputMap",
+        storageClass="HealSparseMap",
+        dimensions=("tract", "patch", "skymap", "band"),
+    )
 
     def __init__(self, *, config=None):
         super().__init__(config=config)
@@ -125,6 +132,9 @@ def __init__(self, *, config=None):
         if not config.doNImage:
             self.outputs.remove("nImage")
 
+        if not self.config.doInputMap:
+            self.outputs.remove("inputMap")
+
 
 class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig,
                           pipelineConnections=AssembleCoaddConnections):
@@ -242,6 +252,15 @@ class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig
         dtype=bool,
         default=False,
     )
+    doInputMap = pexConfig.Field(
+        doc="Create a bitwise map of coadd inputs",
+        dtype=bool,
+        default=False,
+    )
+    inputMapper = pexConfig.ConfigurableField(
+        doc="Input map creation subtask.",
+        target=HealSparseInputMapTask,
+    )
 
     def setDefaults(self):
         super().setDefaults()
@@ -407,6 +426,9 @@ def __init__(self, *args, **kwargs):
                                    mask.getMaskPlaneDict().keys())
             del mask
 
+        if self.config.doInputMap:
+            self.makeSubtask("inputMapper")
+
         self.warpType = self.config.warpType
 
     @utils.inheritDoc(pipeBase.PipelineTask)
@@ -773,6 +795,7 @@ def run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
 
            - ``coaddExposure``: coadded exposure (``lsst.afw.image.Exposure``).
            - ``nImage``: exposure count image (``lsst.afw.image.Image``), if requested.
+           - ``inputMap``: bit-wise map of inputs, if requested.
            - ``warpRefList``: input list of refs to the warps (
                               ``lsst.daf.butler.DeferredDatasetHandle`` or
                               ``lsst.daf.persistence.ButlerDataRef``)
@@ -799,6 +822,13 @@ def run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
             nImage = afwImage.ImageU(skyInfo.bbox)
         else:
             nImage = None
+        # If inputMap is requested, create the initial version that can be masked in
+        # assembleSubregion.
+        if self.config.doInputMap:
+            self.inputMapper.build_ccd_input_map(skyInfo.bbox,
+                                                 skyInfo.wcs,
+                                                 coaddExposure.getInfo().getCoaddInputs().ccds)
+
         for subBBox in self._subBBoxIter(skyInfo.bbox, subregionSize):
             try:
                 self.assembleSubregion(coaddExposure, subBBox, tempExpRefList, imageScalerList,
@@ -807,13 +837,20 @@ def run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
             except Exception as e:
                 self.log.fatal("Cannot compute coadd %s: %s", subBBox, e)
 
+        # If inputMap is requested, we must finalize the map after the accumulation.
+        if self.config.doInputMap:
+            self.inputMapper.finalize_ccd_input_map_mask()
+            inputMap = self.inputMapper.ccd_input_map
+        else:
+            inputMap = None
+
         self.setInexactPsf(coaddMaskedImage.getMask())
         # Despite the name, the following doesn't really deal with "EDGE" pixels: it identifies
         # pixels that didn't receive any unmasked inputs (as occurs around the edge of the field).
         coaddUtils.setCoaddEdgeBits(coaddMaskedImage.getMask(), coaddMaskedImage.getVariance())
         return pipeBase.Struct(coaddExposure=coaddExposure, nImage=nImage,
                                warpRefList=tempExpRefList, imageScalerList=imageScalerList,
-                               weightList=weightList)
+                               weightList=weightList, inputMap=inputMap)
 
     def assembleMetadata(self, coaddExposure, tempExpRefList, weightList):
         """Set the metadata for the coadd.
@@ -947,6 +984,10 @@ def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList
                 self.removeMaskPlanes(maskedImage)
             maskedImageList.append(maskedImage)
 
+            if self.config.doInputMap:
+                visit = exposure.getInfo().getCoaddInputs().visits[0].getId()
+                self.inputMapper.mask_warp_bbox(bbox, visit, mask, statsCtrl.getAndMask())
+
         with self.timer("stack"):
             coaddSubregion = afwMath.statisticsStack(maskedImageList, statsFlags, statsCtrl, weightList,
                                                      clipped,  # also set output to CLIPPED if sigma-clipped
@@ -1545,6 +1586,8 @@ def buildDifferenceImage(self, skyInfo, tempExpRefList, imageScalerList, weightL
         configIntersection = {k: getattr(self.config, k)
                               for k, v in self.config.toDict().items()
                               if (k in config.keys() and k != "connections")}
+        configIntersection['doInputMap'] = False
+        configIntersection['doNImage'] = False
         config.update(**configIntersection)
 
         # statistic MEAN copied from self.config.statistic, but for clarity explicitly assign