Coverage for python/lsst/drp/tasks/forcedPhotCoadd.py: 24%
146 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-08-17 14:28 -0700
1# This file is part of drp_tasks.
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/>.
22import dataclasses
24import lsst.afw.table
25import lsst.pex.config
26import lsst.pipe.base as pipeBase
27from lsst.cell_coadds import MultipleCellCoadd
28from lsst.meas.base._id_generator import SkyMapIdGeneratorConfig
29from lsst.meas.base.applyApCorr import ApplyApCorrTask
30from lsst.meas.base.catalogCalculation import CatalogCalculationTask
31from lsst.meas.base.forcedMeasurement import ForcedMeasurementTask
32from lsst.meas.extensions.scarlet.io import updateCatalogFootprints
34__all__ = ("ForcedPhotCoaddConfig", "ForcedPhotCoaddTask")
37class ForcedPhotCoaddConnections(
38 pipeBase.PipelineTaskConnections,
39 dimensions=("band", "skymap", "tract", "patch"),
40 defaultTemplates={"inputCoaddName": "deep", "outputCoaddName": "deep"},
41):
42 inputSchema = pipeBase.connectionTypes.InitInput(
43 doc="Schema for the input measurement catalogs.",
44 name="{inputCoaddName}Coadd_ref_schema",
45 storageClass="SourceCatalog",
46 )
47 outputSchema = pipeBase.connectionTypes.InitOutput(
48 doc="Schema for the output forced measurement catalogs.",
49 name="{outputCoaddName}Coadd_forced_src_schema",
50 storageClass="SourceCatalog",
51 )
52 exposure = pipeBase.connectionTypes.Input(
53 doc="Input exposure to perform photometry on.",
54 name="{inputCoaddName}Coadd_calexp",
55 storageClass="ExposureF",
56 dimensions=["band", "skymap", "tract", "patch"],
57 )
58 exposure_cell = pipeBase.connectionTypes.Input(
59 doc="Input cell-based coadd exposure to perform photometry on.",
60 name="{inputCoaddName}CoaddCell",
61 storageClass="MultipleCellCoadd",
62 dimensions=["band", "skymap", "tract", "patch"],
63 )
64 background = pipeBase.connectionTypes.Input(
65 doc="Background to subtract from the exposure_cell.",
66 name="{inputCoaddName}Coadd_calexp_background",
67 storageClass="Background",
68 dimensions=["band", "skymap", "tract", "patch"],
69 )
70 refCat = pipeBase.connectionTypes.Input(
71 doc="Catalog of shapes and positions at which to force photometry.",
72 name="{inputCoaddName}Coadd_ref",
73 storageClass="SourceCatalog",
74 dimensions=["skymap", "tract", "patch"],
75 )
76 refCatInBand = pipeBase.connectionTypes.Input(
77 doc="Catalog of shapes and positions in the band having forced photometry done",
78 name="{inputCoaddName}Coadd_meas",
79 storageClass="SourceCatalog",
80 dimensions=("band", "skymap", "tract", "patch"),
81 )
82 footprintCatInBand = pipeBase.connectionTypes.Input(
83 doc="Catalog of footprints to attach to sources",
84 name="{inputCoaddName}Coadd_deblendedFlux",
85 storageClass="SourceCatalog",
86 dimensions=("band", "skymap", "tract", "patch"),
87 )
88 scarletModels = pipeBase.connectionTypes.Input(
89 doc="Multiband scarlet models produced by the deblender",
90 name="{inputCoaddName}Coadd_scarletModelData",
91 storageClass="LsstScarletModelData",
92 dimensions=("tract", "patch", "skymap"),
93 )
94 refWcs = pipeBase.connectionTypes.Input(
95 doc="Reference world coordinate system.",
96 name="{inputCoaddName}Coadd.wcs",
97 storageClass="Wcs",
98 dimensions=["band", "skymap", "tract", "patch"],
99 ) # used in place of a skymap wcs because of DM-28880
100 measCat = pipeBase.connectionTypes.Output(
101 doc="Output forced photometry catalog.",
102 name="{outputCoaddName}Coadd_forced_src",
103 storageClass="SourceCatalog",
104 dimensions=["band", "skymap", "tract", "patch"],
105 )
107 def __init__(self, *, config=None):
108 super().__init__(config=config)
109 if config is None:
110 return
112 if config.footprintDatasetName != "LsstScarletModelData":
113 self.inputs.remove("scarletModels")
114 if config.footprintDatasetName != "DeblendedFlux":
115 self.inputs.remove("footprintCatInBand")
116 if self.config.imageType == "future":
117 self.exposure = dataclasses.replace(self.exposure, storageClass="CellCoadd")
118 del self.exposure_cell
119 del self.background
120 del self.refWcs
121 elif self.config.useCellCoadds:
122 del self.exposure
123 else:
124 del self.exposure_cell
125 del self.background
128class ForcedPhotCoaddConfig(pipeBase.PipelineTaskConfig, pipelineConnections=ForcedPhotCoaddConnections):
129 measurement = lsst.pex.config.ConfigurableField(
130 target=ForcedMeasurementTask, doc="subtask to do forced measurement"
131 )
132 coaddName = lsst.pex.config.Field(
133 doc="coadd name: typically one of deep or goodSeeing",
134 dtype=str,
135 default="deep",
136 )
137 useCellCoadds = lsst.pex.config.Field(
138 doc="Use cell-based coadds for forced measurements?",
139 dtype=bool,
140 default=False,
141 )
142 doApCorr = lsst.pex.config.Field(
143 dtype=bool, default=True, doc="Run subtask to apply aperture corrections"
144 )
145 applyApCorr = lsst.pex.config.ConfigurableField(
146 target=ApplyApCorrTask, doc="Subtask to apply aperture corrections"
147 )
148 catalogCalculation = lsst.pex.config.ConfigurableField(
149 target=CatalogCalculationTask, doc="Subtask to run catalogCalculation plugins on catalog"
150 )
151 footprintDatasetName = lsst.pex.config.Field(
152 doc="Dataset (without coadd prefix) that should be used to obtain (Heavy)Footprints for sources. "
153 "Must have IDs that match those of the reference catalog."
154 "If None, Footprints will be generated by transforming the reference Footprints.",
155 dtype=str,
156 default="LsstScarletModelData",
157 optional=True,
158 )
159 doConserveFlux = lsst.pex.config.Field(
160 dtype=bool,
161 default=True,
162 doc="Whether to use the deblender models as templates to re-distribute the flux "
163 "from the 'exposure' (True), or to perform measurements on the deblender model footprints. "
164 "If footprintDatasetName != 'LsstScarletModelData' then this field is ignored.",
165 )
166 doStripFootprints = lsst.pex.config.Field(
167 dtype=bool,
168 default=True,
169 doc="Whether to strip footprints from the output catalog before "
170 "saving to disk. "
171 "This is usually done when using scarlet models to save disk space.",
172 )
173 hasFakes = lsst.pex.config.Field(
174 dtype=bool,
175 default=False,
176 doc="Should be set to True if fake sources have been inserted into the input data.",
177 )
178 idGenerator = SkyMapIdGeneratorConfig.make_field()
179 imageType = lsst.pex.config.ChoiceField(
180 "Which image type to expect for the input coadd. "
181 "This option only directly affects connection storage classes and hence 'runQuantum'; the 'run' "
182 "method behavior is determined by which type is actually passed in.",
183 allowed={
184 "legacy": (
185 "Read a lsst.cell_coadds.MultipleCellCoadd via 'exposure_cell` and restore 'background' "
186 "(if useCellCoadd) or lsst.afw.image.Exposure via `exposure` (if not useCellCoadd)."
187 ),
188 "future": (
189 "Read lsst.images.cells.CellCoadd via 'exposure'; useCellCoadd is ignored. This choice also "
190 "deletes the 'refWcs' connection, which is always superfluous in all practical usage of this "
191 "task and would need to be set to a different dataset type whenever this choice is used."
192 ),
193 },
194 dtype=str,
195 optional=False,
196 default="legacy",
197 )
199 def setDefaults(self):
200 # Docstring inherited.
201 # Make catalogCalculation a no-op by default as no modelFlux is setup
202 # by default in ForcedMeasurementTask
203 super().setDefaults()
205 self.catalogCalculation.plugins.names = []
206 self.measurement.copyColumns["id"] = "id"
207 self.measurement.copyColumns["parent"] = "parent"
208 self.measurement.plugins.names |= ["base_InputCount", "base_Variance"]
209 self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [
210 "CLIPPED",
211 "SENSOR_EDGE",
212 "REJECTED",
213 "INEXACT_PSF",
214 # TODO DM-44658 and DM-45980: don't have STREAK propagated yet.
215 # "STREAK",
216 ]
217 self.measurement.plugins["base_PixelFlags"].masksFpCenter = [
218 "CLIPPED",
219 "SENSOR_EDGE",
220 "REJECTED",
221 "INEXACT_PSF",
222 # "STREAK",
223 ]
226class ForcedPhotCoaddTask(pipeBase.PipelineTask):
227 """A pipeline task for performing forced measurement on coadd images.
229 Parameters
230 ----------
231 refSchema : `lsst.afw.table.Schema`, optional
232 The schema of the reference catalog, passed to the constructor of the
233 references subtask. Optional, but must be specified if ``initInputs``
234 is not; if both are specified, ``initInputs`` takes precedence.
235 initInputs : `dict`
236 Dictionary that can contain a key ``inputSchema`` containing the
237 schema. If present will override the value of ``refSchema``.
238 **kwds
239 Keyword arguments are passed to the supertask constructor.
240 """
242 ConfigClass = ForcedPhotCoaddConfig
243 _DefaultName = "forcedPhotCoadd"
244 dataPrefix = "deepCoadd_"
246 def __init__(self, refSchema=None, initInputs=None, **kwds):
247 super().__init__(**kwds)
249 if initInputs is not None:
250 refSchema = initInputs["inputSchema"].schema
252 if refSchema is None:
253 raise ValueError("No reference schema provided.")
254 self.makeSubtask("measurement", refSchema=refSchema)
255 # It is necessary to get the schema internal to the forced measurement
256 # task until such a time that the schema is not owned by the
257 # measurement task, but is passed in by an external caller.
258 if self.config.doApCorr:
259 self.makeSubtask("applyApCorr", schema=self.measurement.schema)
260 self.makeSubtask("catalogCalculation", schema=self.measurement.schema)
261 self.outputSchema = lsst.afw.table.SourceCatalog(self.measurement.schema)
263 def runQuantum(self, butlerQC, inputRefs, outputRefs):
264 inputs = butlerQC.get(inputRefs)
266 refCatInBand = inputs.pop("refCatInBand")
267 if self.config.footprintDatasetName == "LsstScarletModelData":
268 footprintData = inputs.pop("scarletModels")
269 elif self.config.footprintDatasetName == "DeblendedFlux":
270 footprintData = inputs.pop("footprintCatIndBand")
271 else:
272 footprintData = None
274 refCat = inputs.pop("refCat")
275 refWcs = inputs.pop("refWcs", None)
277 if self.config.imageType == "future":
278 coadd = inputs.pop("exposure")
279 dataId = inputRefs.exposure.dataId
280 # Instead of going directly from lsst.images.cells.CellCoadd to
281 # Exposure, it's cleaner for now to go through MultipleCellCoadd
282 # because the apCorrMap and ccdInputs need special handling - the
283 # cell-based versions can't be attached to Exposure. Eventually
284 # we'll rewrite the lower-level code to use the lsst.images
285 # equivalents natively.
286 coadd = coadd.to_legacy_cell_coadd()
287 elif self.config.useCellCoadds:
288 coadd = inputs.pop("exposure_cell")
289 dataId = inputRefs.exposure_cell.dataId
290 else:
291 coadd = inputs.pop("exposure")
292 dataId = inputRefs.exposure.dataId
293 if isinstance(coadd, MultipleCellCoadd):
294 stitched_coadd = coadd.stitch()
295 exposure = stitched_coadd.asExposure()
296 if self.config.imageType == "legacy":
297 background = inputs.pop("background")
298 exposure.image -= background.getImage()
299 apCorrMap = stitched_coadd.ap_corr_map
300 else:
301 exposure = coadd
302 apCorrMap = exposure.getInfo().getApCorrMap()
303 if refWcs is None:
304 refWcs = exposure.getWcs()
306 assert not inputs, "runQuantum got extra inputs."
308 measCat, exposureId = self.generateMeasCat(
309 dataId=dataId,
310 exposure=exposure,
311 refCat=refCat,
312 refCatInBand=refCatInBand,
313 refWcs=refWcs,
314 footprintData=footprintData,
315 )
316 outputs = self.run(
317 measCat=measCat,
318 exposure=exposure,
319 refCat=refCat,
320 refWcs=refWcs,
321 exposureId=exposureId,
322 apCorrMap=apCorrMap,
323 )
324 # Strip HeavyFootprints to save space on disk
325 if self.config.footprintDatasetName == "LsstScarletModelData" and self.config.doStripFootprints:
326 sources = outputs.measCat
327 for source in sources[sources["parent"] != 0]:
328 source.setFootprint(None)
329 butlerQC.put(outputs, outputRefs)
331 def generateMeasCat(self, dataId, exposure, refCat, refCatInBand, refWcs, footprintData):
332 """Generate a measurement catalog.
334 Parameters
335 ----------
336 dataId : `lsst.daf.butler.DataCoordinate`
337 Butler data ID for this image, with ``{tract, patch, band}`` keys.
338 exposure : `lsst.afw.image.exposure.Exposure`
339 Exposure to generate the catalog for.
340 refCat : `lsst.afw.table.SourceCatalog`
341 Catalog of shapes and positions at which to force photometry.
342 refCatInBand : `lsst.afw.table.SourceCatalog`
343 Catalog of shapes and position in the band forced photometry is
344 currently being performed
345 refWcs : `lsst.afw.image.SkyWcs`
346 Reference world coordinate system.
347 footprintData : `ScarletDataModel` or `lsst.afw.table.SourceCatalog`
348 Either the scarlet data models or the deblended catalog containings
349 footprints. If `footprintData` is `None` then the footprints
350 contained in `refCatInBand` are used.
352 Returns
353 -------
354 measCat : `lsst.afw.table.SourceCatalog`
355 Catalog of forced sources to measure.
356 expId : `int`
357 Unique binary id associated with the input exposure
359 Raises
360 ------
361 LookupError
362 Raised if a footprint with a given source id was in the reference
363 catalog but not in the reference catalog in band (meaning there was
364 some sort of mismatch in the two input catalogs)
365 """
366 id_generator = self.config.idGenerator.apply(dataId)
367 measCat = self.measurement.generateMeasCat(
368 exposure, refCat, refWcs, idFactory=id_generator.make_table_id_factory()
369 )
370 # attach footprints here as this can naturally live inside this method
371 if self.config.footprintDatasetName == "LsstScarletModelData":
372 # Load the scarlet models
373 self._attachScarletFootprints(
374 catalog=measCat, modelData=footprintData, exposure=exposure, band=dataId["band"]
375 )
376 else:
377 if self.config.footprintDatasetName is None:
378 footprintCat = refCatInBand
379 else:
380 footprintCat = footprintData
381 for srcRecord in measCat:
382 fpRecord = footprintCat.find(srcRecord.getId())
383 if fpRecord is None:
384 raise LookupError(
385 "Cannot find Footprint for source {}; please check that {} "
386 "IDs are compatible with reference source IDs".format(srcRecord.getId(), footprintCat)
387 )
388 srcRecord.setFootprint(fpRecord.getFootprint())
389 return measCat, id_generator.catalog_id
391 def run(self, measCat, exposure, refCat, refWcs, exposureId=None, apCorrMap=None):
392 """Perform forced measurement on a single exposure.
394 Parameters
395 ----------
396 measCat : `lsst.afw.table.SourceCatalog`
397 The measurement catalog, based on the sources listed in the
398 reference catalog.
399 exposure : `lsst.afw.image.Exposure`
400 The measurement image upon which to perform forced detection.
401 refCat : `lsst.afw.table.SourceCatalog`
402 The reference catalog of sources to measure.
403 refWcs : `lsst.afw.image.SkyWcs`
404 The WCS for the references.
405 exposureId : `int`
406 Optional unique exposureId used for random seed in measurement
407 task.
408 apCorrMap : `~lsst.afw.image.ApCorrMap`, optional
409 Aperture correction map to use for aperture corrections.
410 If not provided, the map is read from the exposure.
412 Returns
413 -------
414 result : ~`lsst.pipe.base.Struct`
415 Structure with fields:
417 ``measCat``
418 Catalog of forced measurement results
419 (`lsst.afw.table.SourceCatalog`).
420 """
421 # We want to cache repeated PSF evaluations at the same point coming
422 # from different measurement plugins. We assume each algorithm tries
423 # to evaluate the PSF twice, which is more than enough since many don't
424 # evaluate it at all, and there's no *good* reason for any algorithm to
425 # evaluate it more than once.
426 exposure.psf.setCacheCapacity(2 * len(self.config.measurement.plugins.names))
427 # Some mask planes may not be defined on the coadds always.
428 # We add the mask planes, which is a no-op if already defined.
429 for maskPlane in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
430 exposure.mask.addMaskPlane(maskPlane)
431 for maskPlane in self.config.measurement.plugins["base_PixelFlags"].masksFpCenter:
432 exposure.mask.addMaskPlane(maskPlane)
433 self.measurement.run(measCat, exposure, refCat, refWcs, exposureId=exposureId)
434 if self.config.doApCorr:
435 if apCorrMap is None:
436 apCorrMap = exposure.getInfo().getApCorrMap()
437 self.applyApCorr.run(catalog=measCat, apCorrMap=apCorrMap)
439 self.catalogCalculation.run(measCat)
441 return pipeBase.Struct(measCat=measCat)
443 def _attachScarletFootprints(self, catalog, modelData, exposure, band):
444 """Attach scarlet models as HeavyFootprints"""
445 if self.config.doConserveFlux:
446 redistributeImage = exposure
447 else:
448 redistributeImage = None
449 # Attach the footprints
450 updateCatalogFootprints(
451 modelData=modelData,
452 catalog=catalog,
453 band=band,
454 imageForRedistribution=redistributeImage,
455 removeScarletData=True,
456 updateFluxColumns=False,
457 )