From b2f7d33a0a94bcbda52c6881c5bd750428181308 Mon Sep 17 00:00:00 2001 From: James Mullaney Date: Mon, 16 Mar 2026 23:34:53 +0000 Subject: [PATCH 1/2] Add task to determine which PVIs overlap a patch --- python/lsst/analysis/tools/tasks/__init__.py | 1 + .../tools/tasks/coaddInputAnalysisTask.py | 229 ++++++++++++++++++ 2 files changed, 230 insertions(+) create mode 100644 python/lsst/analysis/tools/tasks/coaddInputAnalysisTask.py diff --git a/python/lsst/analysis/tools/tasks/__init__.py b/python/lsst/analysis/tools/tasks/__init__.py index 5d33cd197..31f38cbbc 100644 --- a/python/lsst/analysis/tools/tasks/__init__.py +++ b/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 * diff --git a/python/lsst/analysis/tools/tasks/coaddInputAnalysisTask.py b/python/lsst/analysis/tools/tasks/coaddInputAnalysisTask.py new file mode 100644 index 000000000..3f72fd12d --- /dev/null +++ b/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 . +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) From 9d23c5d227daaab978d3747f204f8099c6c49dff Mon Sep 17 00:00:00 2001 From: jrmullaney Date: Tue, 17 Mar 2026 18:31:15 +0000 Subject: [PATCH 2/2] Add atool to calc fraction of PVIs in coadd --- .../analysis/tools/atools/coaddInputCount.py | 54 ++++++++++++++++++- 1 file changed, 52 insertions(+), 2 deletions(-) diff --git a/python/lsst/analysis/tools/atools/coaddInputCount.py b/python/lsst/analysis/tools/atools/coaddInputCount.py index 9374d598d..130425ff3 100644 --- a/python/lsst/analysis/tools/atools/coaddInputCount.py +++ b/python/lsst/analysis/tools/atools/coaddInputCount.py @@ -20,14 +20,22 @@ # along with this program. If not, see . 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": "", + }