DM-54368: Make a metric that records the numbers of visit-detectors that were rejected from a coadd.

python/lsst/analysis/tools/atools/coaddInputCount.py

@@ -20,14 +20,22 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 from __future__ import annotations
 
-__all__ = ("CoaddInputCount", "CoaddQualityCheck", "CoaddQualityPlot")
+__all__ = ("CoaddInputCount", "CoaddQualityCheck", "CoaddQualityPlot", "CoaddInputFraction")
 
 from lsst.pex.config import ListField
 
 from ..actions.plot.calculateRange import MinMax
 from ..actions.plot.coaddDepthPlot import CoaddDepthPlot
 from ..actions.plot.skyPlot import SkyPlot
-from ..actions.scalar.scalarActions import MeanAction, MedianAction, SigmaMadAction, StdevAction
+from ..actions.scalar.scalarActions import (
+    CountAction,
+    DivideScalar,
+    MeanAction,
+    MedianAction,
+    SigmaMadAction,
+    StdevAction,
+    SumAction,
+)
 from ..actions.vector import BandSelector, CoaddPlotFlagSelector, DownselectVector, LoadVector, SnSelector
 from ..interfaces import AnalysisTool
 
@@ -213,3 +221,45 @@ def setDefaults(self):
         self.process.buildActions.pixels = LoadVector(vectorKey="pixels")
 
         self.produce.plot = CoaddDepthPlot()
+
+
+class CoaddInputFraction(AnalysisTool):
+    """Metrics quantifying the number of images that overlap a patch,
+    the number that actually made it into the coadd, and the ratio
+    of the two (as a fraction).
+    """
+
+    parameterizedBand: bool = False
+
+    def setDefaults(self):
+        super().setDefaults()
+
+        # The number of entries is the number of potential raws:
+        self.process.calculateActions.rawCount = CountAction(vectorKey="inCoadd")
+        self.process.calculateActions.overlapCount = SumAction(vectorKey="patchOverlap")
+        self.process.calculateActions.inVisitSummaryCount = SumAction(vectorKey="visitSummaryRecord")
+        self.process.calculateActions.inCoaddCount = SumAction(vectorKey="inCoadd")
+
+        self.process.calculateActions.inVisitSummaryFraction = DivideScalar()
+        self.process.calculateActions.inVisitSummaryFraction.actionA = SumAction(
+            vectorKey="visitSummaryRecord"
+        )
+        self.process.calculateActions.inVisitSummaryFraction.actionB = CountAction(vectorKey="inCoadd")
+
+        self.process.calculateActions.overlapFraction = DivideScalar()
+        self.process.calculateActions.overlapFraction.actionA = SumAction(vectorKey="patchOverlap")
+        self.process.calculateActions.overlapFraction.actionB = CountAction(vectorKey="inCoadd")
+
+        self.process.calculateActions.inCoaddFraction = DivideScalar()
+        self.process.calculateActions.inCoaddFraction.actionA = SumAction(vectorKey="inCoadd")
+        self.process.calculateActions.inCoaddFraction.actionB = CountAction(vectorKey="inCoadd")
+
+        self.produce.metric.units = {
+            "rawCount": "ct",
+            "inVisitSummaryCount": "ct",
+            "overlapCount": "ct",
+            "inCoaddCount": "ct",
+            "inVisitSummaryFraction": "",
+            "overlapFraction": "",
+            "inCoaddFraction": "",
+        }

python/lsst/analysis/tools/tasks/__init__.py

@@ -9,6 +9,7 @@
 from .coaddDepthSummary import *
 from .coaddDepthSummaryPlot import *
 from .coaddDepthTableTractAnalysis import *
+from .coaddInputAnalysisTask import *
 from .diaFakesDetectorVisitAnalysis import *
 from .diaFakesVisitAnalysis import *
 from .diaObjectDetectorVisitAnalysis import *

python/lsst/analysis/tools/tasks/coaddInputAnalysisTask.py

@@ -0,0 +1,229 @@
+# This file is part of analysis_tools.
+#
+# 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/>.
+from __future__ import annotations
+
+__all__ = (
+    "CoaddInputAnalysisConfig",
+    "CoaddInputAnalysisTask",
+)
+
+from collections import defaultdict
+
+import astropy.table
+import numpy as np
+
+import lsst.geom
+import lsst.sphgeom as sphgeom
+from lsst.pipe.base import (
+    InputQuantizedConnection,
+    OutputQuantizedConnection,
+    QuantumContext,
+)
+from lsst.pipe.base import connectionTypes as cT
+from lsst.skymap import BaseSkyMap
+
+from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask
+
+
+class CoaddInputAnalysisConnections(
+    AnalysisBaseConnections,
+    dimensions=("tract", "patch", "band"),
+    defaultTemplates={"inputName": "raw", "outputName": "coaddInputAnalysis"},
+):
+    raw = cT.Input(
+        doc="Raw exposure data.",
+        name="{inputName}",
+        storageClass="Exposure",
+        dimensions=(
+            "instrument",
+            "exposure",
+            "detector",
+        ),
+        deferLoad=True,
+        multiple=True,
+    )
+    coaddInputs = cT.Input(
+        doc="List of dataIds that went into a coadd.",
+        name="deepCoadd_input_summary_tract",
+        storageClass="ArrowAstropy",
+        dimensions=(
+            "skymap",
+            "tract",
+            "band",
+        ),
+    )
+    visitTable = cT.Input(
+        doc="""Table summarising the visit images.
+            "Contains image corners that are used to
+            "refine the coverage.""",
+        name="visit_detector_table",
+        storageClass="ArrowAstropy",
+        dimensions=("instrument",),
+        deferLoad=True,
+        multiple=False,
+    )
+    # visitSummary isn't used, but a per-visit data input
+    # substantially reduces the time to build quantum graph.
+    visitSummary = cT.Input(
+        doc="""Table summarising the visit images.
+            "Contains image corners that are used to
+            "refine the coverage.""",
+        name="visit_summary",
+        storageClass="ExposureCatalog",
+        dimensions=(
+            "instrument",
+            "visit",
+        ),
+        deferLoad=True,
+        multiple=True,
+    )
+    skyMap = cT.PrerequisiteInput(
+        doc="Sky map defining the tracts and patches.",
+        name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
+        storageClass="SkyMap",
+        dimensions=("skymap",),
+    )
+
+
+class CoaddInputAnalysisConfig(AnalysisBaseConfig, pipelineConnections=CoaddInputAnalysisConnections):
+    pass
+
+
+class CoaddInputAnalysisTask(AnalysisPipelineTask):
+    """
+    Construct a table containing the visit and detector IDs
+    of all PVIs that the butler believes could potentially
+    have gone into a coadd. Indicate those whose post-calibration
+    astrometry overlaps with the patch. Also Indicate those which
+    finally made it into the coadd.
+    """
+
+    ConfigClass = CoaddInputAnalysisConfig
+    _DefaultName = "coaddInputAnalysis"
+
+    def runQuantum(
+        self,
+        butlerQC: QuantumContext,
+        inputRefs: InputQuantizedConnection,
+        outputRefs: OutputQuantizedConnection,
+    ) -> None:
+
+        # The PVIs that actually made it into the coadd:
+        coaddInputTable = butlerQC.get(inputRefs.coaddInputs)
+        patch = butlerQC.quantum.dataId["patch"]
+        patchInputs = coaddInputTable[coaddInputTable["patch"] == patch]
+        inCoadd = set(zip(patchInputs["visit"], patchInputs["detector"]))
+
+        # On-sky polygon of patch to determine which PVIs cover it.
+        dataId = butlerQC.quantum.dataId
+        skyMap = butlerQC.get(inputRefs.skyMap)
+        tractInfo = skyMap[dataId["tract"]]
+        patchInfo = tractInfo.getPatchInfo(dataId["patch"])
+        patchPoly = patchInfo.getOuterSkyPolygon()
+
+        visitTableHandle = butlerQC.get(inputRefs.visitTable)
+        columns = [
+            "visitId",
+            "detector",
+            "llcra",
+            "llcdec",
+            "ulcra",
+            "ulcdec",
+            "urcra",
+            "urcdec",
+            "lrcra",
+            "lrcdec",
+        ]
+        imageCorners = self.loadData(visitTableHandle, columns)
+
+        # Group raws by visit so the imageCorners table can be reduced by visit
+        rawsByVisit = defaultdict(list)
+        for rawRef in inputRefs.raw:
+            rawsByVisit[rawRef.dataId["exposure"]].append(rawRef)
+
+        visits = []
+        detectors = []
+        inCoadd_col = []
+        visitSummaryRecord = []
+        overlapsWithPatch = []
+        for visit, rawRefs in rawsByVisit.items():
+            visitSummary = imageCorners[imageCorners["visitId"] == visit]
+
+            if len(visitSummary) == 0:
+                visits += [visit] * len(rawRefs)
+                detectors += [rawRef.dataId["detector"] for rawRef in rawRefs]
+                inCoadd_col += [False] * len(rawRefs)
+                visitSummaryRecord += [False] * len(rawRefs)
+                overlapsWithPatch += [False] * len(rawRefs)
+                continue
+
+            for rawRef in rawRefs:
+                detector = rawRef.dataId["detector"]
+                visits.append(visit)
+                detectors.append(detector)
+                inCoadd_col.append((visit, detector) in inCoadd)
+
+                if not np.any(rowMask := visitSummary["detector"] == detector):
+                    # If there's no visit summary, record and move on:
+                    visitSummaryRecord.append(False)
+                    overlapsWithPatch.append(False)
+                    continue
+
+                # If we have got this far, then there is a row for this image.
+                # Record whether the calibrated image covers the patch:
+                visitSummaryRecord.append(True)
+
+                raCorners = np.array(
+                    [
+                        visitSummary[rowMask][corner].value[0]
+                        for corner in ["llcra", "ulcra", "urcra", "lrcra"]
+                    ]
+                )
+                decCorners = np.array(
+                    [
+                        visitSummary[rowMask][corner].value[0]
+                        for corner in ["llcdec", "ulcdec", "urcdec", "lrcdec"]
+                    ]
+                )
+
+                if np.all(np.isfinite(raCorners)) and np.all(np.isfinite(decCorners)):
+                    detCorners = [
+                        lsst.geom.SpherePoint(ra, dec, units=lsst.geom.degrees).getVector()
+                        for ra, dec in zip(raCorners, decCorners)
+                    ]
+                    detPoly = sphgeom.ConvexPolygon.convexHull(detCorners)
+                    overlapsWithPatch.append(patchPoly.intersects(detPoly))
+                else:
+                    overlapsWithPatch.append(False)
+
+        data = astropy.table.Table(
+            {
+                "visit": visits,
+                "detector": detectors,
+                "visitSummaryRecord": visitSummaryRecord,
+                "patchOverlap": overlapsWithPatch,
+                "inCoadd": inCoadd_col,
+            }
+        )
+
+        # Send the data table to the atools for analysis:
+        outputs = self.run(data=data)
+        butlerQC.put(outputs, outputRefs)