Coverage for python / lsst / analysis / tools / tasks / coaddInputAnalysis.py: 26%
76 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-06 09:07 +0000
1# This file is part of analysis_tools.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21from __future__ import annotations
23__all__ = (
24 "CoaddInputAnalysisConfig",
25 "CoaddInputAnalysisTask",
26)
28from collections import defaultdict
30import astropy.table
31import numpy as np
33from lsst.geom import SpherePoint, degrees
34from lsst.pipe.base import (
35 InputQuantizedConnection,
36 OutputQuantizedConnection,
37 QuantumContext,
38)
39from lsst.pipe.base import connectionTypes as cT
40from lsst.skymap import BaseSkyMap
41from lsst.sphgeom import ConvexPolygon
43from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask
46class CoaddInputAnalysisConnections(
47 AnalysisBaseConnections,
48 dimensions=("tract", "patch", "band"),
49 defaultTemplates={"outputName": "coaddInputAnalysis"},
50):
51 raw = cT.Input(
52 doc="Raw exposure data.",
53 name="raw",
54 storageClass="Exposure",
55 dimensions=(
56 "instrument",
57 "exposure",
58 "detector",
59 ),
60 deferLoad=True,
61 multiple=True,
62 )
63 coaddInputs = cT.Input(
64 doc="List of dataIds that went into a coadd.",
65 name="deep_coadd_input_summary_tract",
66 storageClass="ArrowAstropy",
67 dimensions=(
68 "skymap",
69 "tract",
70 "band",
71 ),
72 )
73 visitTable = cT.Input(
74 doc="""Table summarising metadata at detector-level for all visist,
75 from which the per-detector image corner coordinates are used
76 to determine which calibrated images spatially overlap a given
77 patch.""",
78 name="ccdVisitTable",
79 storageClass="ArrowAstropy",
80 dimensions=("instrument",),
81 deferLoad=True,
82 multiple=False,
83 )
84 visitSummary = cT.Input(
85 doc="""Per-visit summary of calibrated images. Not directly used in
86 processing, but including it as a per-visit input substantially
87 reduces the time to build the quantum graph.""",
88 name="finalVisitSummary",
89 storageClass="ExposureCatalog",
90 dimensions=(
91 "instrument",
92 "visit",
93 ),
94 deferLoad=True,
95 multiple=True,
96 )
97 skyMap = cT.PrerequisiteInput(
98 doc="Sky map defining the tracts and patches.",
99 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
100 storageClass="SkyMap",
101 dimensions=("skymap",),
102 )
105class CoaddInputAnalysisConfig(AnalysisBaseConfig, pipelineConnections=CoaddInputAnalysisConnections):
106 pass
109class CoaddInputAnalysisTask(AnalysisPipelineTask):
110 """
111 Construct a table containing the visit and detector IDs
112 of all PVIs that the butler believes could potentially
113 have gone into a coadd. Indicate those whose post-calibration
114 astrometry overlaps with the patch. Also Indicate those which
115 finally made it into the coadd.
116 """
118 ConfigClass = CoaddInputAnalysisConfig
119 _DefaultName = "coaddInputAnalysis"
121 def runQuantum(
122 self,
123 butlerQC: QuantumContext,
124 inputRefs: InputQuantizedConnection,
125 outputRefs: OutputQuantizedConnection,
126 ) -> None:
128 # The PVIs that actually made it into the coadd:
129 coaddInputTable = butlerQC.get(inputRefs.coaddInputs)
130 patch = butlerQC.quantum.dataId["patch"]
131 patchInputs = coaddInputTable[coaddInputTable["patch"] == patch]
132 inCoadd = set(zip(patchInputs["visit"], patchInputs["detector"]))
134 # On-sky polygon of patch to determine which PVIs cover it.
135 dataId = butlerQC.quantum.dataId
136 skyMap = butlerQC.get(inputRefs.skyMap)
137 tractInfo = skyMap[dataId["tract"]]
138 patchInfo = tractInfo.getPatchInfo(dataId["patch"])
139 patchPoly = patchInfo.getOuterSkyPolygon()
141 visitTableHandle = butlerQC.get(inputRefs.visitTable)
142 columns = [
143 "visitId",
144 "detector",
145 "llcra",
146 "llcdec",
147 "ulcra",
148 "ulcdec",
149 "urcra",
150 "urcdec",
151 "lrcra",
152 "lrcdec",
153 ]
154 imageCorners = self.loadData(visitTableHandle, columns)
156 # Group raws by visit so the imageCorners table can be reduced by visit
157 rawsByVisit = defaultdict(list)
158 for rawRef in inputRefs.raw:
159 rawsByVisit[rawRef.dataId["exposure"]].append(rawRef)
161 data = self.makeData(inCoadd, patchPoly, imageCorners, rawsByVisit)
163 outputs = self.run(data=data)
164 butlerQC.put(outputs, outputRefs)
166 def makeData(self, inCoadd, patchPoly, imageCorners, rawsByVisit):
167 """Build the per-PVI data table for analysis.
169 Parameters
170 ----------
171 inCoadd : `set` [`tuple` [`int`, `int`]]
172 Set of (visit, detector) pairs that made it into the coadd.
173 patchPoly : `lsst.sphgeom.ConvexPolygon`
174 On-sky polygon of the patch.
175 imageCorners : `astropy.table.Table`
176 Table of per-detector image corner coordinates.
177 rawsByVisit : `dict` [`int`, `list`]
178 Raw data references grouped by visit ID.
180 Returns
181 -------
182 data : `astropy.table.Table`
183 Table with columns ``visit``, ``detector``,
184 ``visitSummaryRecord``, ``patchOverlap``, and ``inCoadd``.
185 """
186 n = sum(len(rawRefs) for rawRefs in rawsByVisit.values())
187 visits = [None] * n
188 detectors = [None] * n
189 inCoadd_col = [False] * n
190 visitSummaryRecord = [False] * n
191 overlapsWithPatch = [False] * n
193 i = 0
194 for visit, rawRefs in rawsByVisit.items():
195 visitSummary = imageCorners[imageCorners["visitId"] == visit]
197 if len(visitSummary) == 0:
198 for rawRef in rawRefs:
199 visits[i] = visit
200 detectors[i] = rawRef.dataId["detector"]
201 i += 1
202 continue
204 for rawRef in rawRefs:
205 detector = rawRef.dataId["detector"]
206 visits[i] = visit
207 detectors[i] = detector
208 inCoadd_col[i] = (visit, detector) in inCoadd
210 if not np.any(rowMask := visitSummary["detector"] == detector):
211 # If there's no visit summary, move on.
212 # visitSummaryRecord and overlapsWithPatch already contain
213 # False as default:
214 i += 1
215 continue
217 # If we have got this far, then there is a row for this image.
218 # Record whether the calibrated image covers the patch:
219 visitSummaryRecord[i] = True
221 # It is a Bad Thing if there is more than one row per
222 # (visit,detector):
223 if (nMatches := rowMask.sum()) > 1:
224 raise RuntimeError(
225 f"Expected exactly one row for visit={visit}, "
226 f"detector={detector} in imageCorners, found {nMatches}."
227 )
228 raCorners = np.array(
229 [
230 visitSummary[rowMask][corner].value[0]
231 for corner in ["llcra", "ulcra", "urcra", "lrcra"]
232 ]
233 )
234 decCorners = np.array(
235 [
236 visitSummary[rowMask][corner].value[0]
237 for corner in ["llcdec", "ulcdec", "urcdec", "lrcdec"]
238 ]
239 )
241 # Defensively catch non-finite corner coords. Treat such
242 # detectors conservatively as not overlapping the patch.
243 if np.all(np.isfinite(raCorners)) and np.all(np.isfinite(decCorners)):
244 detCorners = [
245 SpherePoint(ra, dec, units=degrees).getVector()
246 for ra, dec in zip(raCorners, decCorners)
247 ]
248 detPoly = ConvexPolygon.convexHull(detCorners)
249 overlapsWithPatch[i] = patchPoly.intersects(detPoly)
251 i += 1
253 return astropy.table.Table(
254 {
255 "visit": visits,
256 "detector": detectors,
257 "visitSummaryRecord": visitSummaryRecord,
258 "patchOverlap": overlapsWithPatch,
259 "inCoadd": inCoadd_col,
260 }
261 )