Coverage for python/lsst/pipe/tasks/multiBand.py: 22%
349 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-19 07:35 +0000
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-19 07:35 +0000
1# This file is part of pipe_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/>.
22__all__ = ["DetectCoaddSourcesConfig", "DetectCoaddSourcesTask",
23 "MeasureMergedCoaddSourcesConfig", "MeasureMergedCoaddSourcesTask",
24 ]
26import numpy as np
28from lsst.geom import Extent2I
29from lsst.pipe.base import (
30 AnnotatedPartialOutputsError,
31 Struct,
32 PipelineTask,
33 PipelineTaskConfig,
34 PipelineTaskConnections
35)
36import lsst.pipe.base.connectionTypes as cT
37from lsst.pex.config import Field, ConfigurableField, ChoiceField
38from lsst.meas.algorithms import (
39 DynamicDetectionTask,
40 ExceedsMaxVarianceScaleError,
41 InsufficientSourcesError,
42 PsfGenerationError,
43 ReferenceObjectLoader,
44 ScaleVarianceTask,
45 SetPrimaryFlagsTask,
46 TooManyMaskedPixelsError,
47 ZeroFootprintError,
48)
49from lsst.meas.base import (
50 SingleFrameMeasurementTask,
51 ApplyApCorrTask,
52 CatalogCalculationTask,
53 SkyMapIdGeneratorConfig,
54)
55from lsst.meas.extensions.scarlet.io import updateCatalogFootprints
56from lsst.meas.astrom import DirectMatchTask, denormalizeMatches
57from lsst.pipe.tasks.propagateSourceFlags import PropagateSourceFlagsTask
58import lsst.afw.image as afwImage
59import lsst.afw.math as afwMath
60import lsst.afw.table as afwTable
61from lsst.daf.base import PropertyList
62from lsst.skymap import BaseSkyMap
64# NOTE: these imports are a convenience so multiband users only have to import this file.
65from .mergeDetections import MergeDetectionsConfig, MergeDetectionsTask # noqa: F401
66from .mergeMeasurements import MergeMeasurementsConfig, MergeMeasurementsTask # noqa: F401
67from .multiBandUtils import CullPeaksConfig # noqa: F401
68from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiConfig # noqa: F401
69from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiTask # noqa: F401
72"""
73New set types:
74* deepCoadd_det: detections from what used to be processCoadd (tract, patch, filter)
75* deepCoadd_mergeDet: merged detections (tract, patch)
76* deepCoadd_meas: measurements of merged detections (tract, patch, filter)
77* deepCoadd_ref: reference sources (tract, patch)
78All of these have associated *_schema catalogs that require no data ID and hold no records.
80In addition, we have a schema-only dataset, which saves the schema for the PeakRecords in
81the mergeDet, meas, and ref dataset Footprints:
82* deepCoadd_peak_schema
83"""
86##############################################################################################################
87class DetectCoaddSourcesConnections(PipelineTaskConnections,
88 dimensions=("tract", "patch", "band", "skymap"),
89 defaultTemplates={"inputCoaddName": "deep", "outputCoaddName": "deep"}):
90 detectionSchema = cT.InitOutput(
91 doc="Schema of the detection catalog",
92 name="{outputCoaddName}Coadd_det_schema",
93 storageClass="SourceCatalog",
94 )
95 exposure = cT.Input(
96 doc="Exposure on which detections are to be performed. ",
97 name="{inputCoaddName}Coadd",
98 storageClass="ExposureF",
99 dimensions=("tract", "patch", "band", "skymap")
100 )
101 exposure_cells = cT.Input(
102 doc="Exposure on which detections are to be performed. ",
103 name="{inputCoaddName}CoaddCell",
104 storageClass="MultipleCellCoadd",
105 dimensions=("tract", "patch", "band", "skymap"),
106 )
107 skyMap = cT.Input(
108 doc="Description of the skymap's tracts and patches.",
109 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
110 storageClass="SkyMap",
111 dimensions=("skymap",),
112 )
113 outputBackgrounds = cT.Output(
114 doc="Output Backgrounds used in detection",
115 name="{outputCoaddName}Coadd_calexp_background",
116 storageClass="Background",
117 dimensions=("tract", "patch", "band", "skymap")
118 )
119 outputSources = cT.Output(
120 doc="Detected sources catalog",
121 name="{outputCoaddName}Coadd_det",
122 storageClass="SourceCatalog",
123 dimensions=("tract", "patch", "band", "skymap")
124 )
125 outputExposure = cT.Output(
126 doc="Exposure post detection",
127 name="{outputCoaddName}Coadd_calexp",
128 storageClass="ExposureF",
129 dimensions=("tract", "patch", "band", "skymap")
130 )
132 def __init__(self, *, config=None):
133 super().__init__(config=config)
134 assert isinstance(config, DetectCoaddSourcesConfig)
136 if config.useCellCoadds:
137 del self.exposure
138 else:
139 del self.exposure_cells
141 if not self.config.forceExactBinning:
142 del self.skyMap
143 if self.config.writeOnlyBackgrounds:
144 del self.outputExposure
145 del self.outputSources
146 del self.detectionSchema
149class DetectCoaddSourcesConfig(PipelineTaskConfig, pipelineConnections=DetectCoaddSourcesConnections):
150 """Configuration parameters for the DetectCoaddSourcesTask
151 """
153 doScaleVariance = Field(dtype=bool, default=True, doc="Scale variance plane using empirical noise?")
154 scaleVariance = ConfigurableField(target=ScaleVarianceTask, doc="Variance rescaling")
155 detection = ConfigurableField(target=DynamicDetectionTask, doc="Source detection")
156 coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
157 useCellCoadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?")
158 hasFakes = Field(
159 dtype=bool,
160 default=False,
161 doc="Should be set to True if fake sources have been inserted into the input data.",
162 )
163 idGenerator = SkyMapIdGeneratorConfig.make_field()
164 forceExactBinning = Field(
165 dtype=bool,
166 default=False,
167 doc=(
168 "Check that the background bin size evenly divides the patch inner region, and "
169 "crop the outer region to an integer number of bins."
170 )
171 )
172 writeOnlyBackgrounds = Field(dtype=bool, default=False, doc="If true, only save the background models.")
173 writeEmptyBackgrounds = Field(
174 dtype=bool,
175 default=True,
176 doc=(
177 "If true, save a placeholder background with NaNs in all bins (but the right geometry) when "
178 "there are no pixels to compute a background from. This can be useful if a later task combines "
179 "backgrounds from multiple patches as input."
180 )
181 )
183 def setDefaults(self):
184 super().setDefaults()
185 self.detection.thresholdType = "pixel_stdev"
186 self.detection.isotropicGrow = True
187 # Coadds are made from background-subtracted CCDs, so any background subtraction should be very basic
188 self.detection.reEstimateBackground = False
189 self.detection.background.useApprox = False
190 self.detection.background.binSize = 4096
191 self.detection.background.undersampleStyle = 'REDUCE_INTERP_ORDER'
192 self.detection.doTempWideBackground = True # Suppress large footprints that overwhelm the deblender
193 # Include band in packed data IDs that go into object IDs (None -> "as
194 # many bands as are defined", rather than the default of zero).
195 self.idGenerator.packer.n_bands = None
198class DetectCoaddSourcesTask(PipelineTask):
199 """Detect sources on a single filter coadd.
201 Coadding individual visits requires each exposure to be warped. This
202 introduces covariance in the noise properties across pixels. Before
203 detection, we correct the coadd variance by scaling the variance plane in
204 the coadd to match the observed variance. This is an approximate
205 approach -- strictly, we should propagate the full covariance matrix --
206 but it is simple and works well in practice.
208 After scaling the variance plane, we detect sources and generate footprints
209 by delegating to the @ref SourceDetectionTask_ "detection" subtask.
211 DetectCoaddSourcesTask is meant to be run after assembling a coadded image
212 in a given band. The purpose of the task is to update the background,
213 detect all sources in a single band and generate a set of parent
214 footprints. Subsequent tasks in the multi-band processing procedure will
215 merge sources across bands and, eventually, perform forced photometry.
217 Parameters
218 ----------
219 schema : `lsst.afw.table.Schema`, optional
220 Initial schema for the output catalog, modified-in place to include all
221 fields set by this task. If None, the source minimal schema will be used.
222 **kwargs
223 Additional keyword arguments.
224 """
226 _DefaultName = "detectCoaddSources"
227 ConfigClass = DetectCoaddSourcesConfig
229 def __init__(self, schema=None, **kwargs):
230 # N.B. Super is used here to handle the multiple inheritance of PipelineTasks, the init tree
231 # call structure has been reviewed carefully to be sure super will work as intended.
232 super().__init__(**kwargs)
233 if schema is None:
234 schema = afwTable.SourceTable.makeMinimalSchema()
235 self.schema = schema
236 self.makeSubtask("detection", schema=self.schema)
237 if self.config.doScaleVariance:
238 self.makeSubtask("scaleVariance")
240 self.detectionSchema = afwTable.SourceCatalog(self.schema)
242 def runQuantum(self, butlerQC, inputRefs, outputRefs):
243 inputs = butlerQC.get(inputRefs)
244 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
246 if self.config.useCellCoadds:
247 multiple_cell_coadd = inputs.pop("exposure_cells")
248 exposure = multiple_cell_coadd.stitch().asExposure()
249 else:
250 exposure = inputs.pop("exposure")
252 skyMap = inputs.pop("skyMap", None)
253 if skyMap is not None:
254 patchInfo = skyMap[butlerQC.quantum.dataId["tract"]][butlerQC.quantum.dataId["patch"]]
255 else:
256 patchInfo = None
258 assert not inputs, "runQuantum got more inputs than expected."
259 try:
260 outputs = self.run(
261 exposure=exposure,
262 idFactory=idGenerator.make_table_id_factory(),
263 expId=idGenerator.catalog_id,
264 patchInfo=patchInfo,
265 )
266 except (
267 TooManyMaskedPixelsError,
268 ExceedsMaxVarianceScaleError,
269 InsufficientSourcesError,
270 PsfGenerationError,
271 ZeroFootprintError,
272 ) as e:
273 if self.config.writeEmptyBackgrounds:
274 butlerQC.put(self._makeEmptyBackground(exposure, patchInfo), outputRefs.outputBackgrounds)
275 # Detection failed, so clear any leftover the detected mask planes.
276 for maskName in ["DETECTED", "DETECTED_NEGATIVE"]:
277 if maskName in exposure.mask.getMaskPlaneDict().keys():
278 detectedMask = exposure.mask.getMaskPlane(maskName)
279 exposure.mask.clearMaskPlane(detectedMask)
280 butlerQC.put(exposure, outputRefs.outputExposure)
281 error = AnnotatedPartialOutputsError.annotate(
282 e,
283 self,
284 exposure,
285 log=self.log,
286 )
287 raise error from e
289 butlerQC.put(outputs, outputRefs)
291 def run(self, exposure, idFactory, expId, patchInfo=None):
292 """Run detection on an exposure.
294 First scale the variance plane to match the observed variance
295 using ``ScaleVarianceTask``. Then invoke the ``SourceDetectionTask_`` "detection" subtask to
296 detect sources.
298 Parameters
299 ----------
300 exposure : `lsst.afw.image.Exposure`
301 Exposure on which to detect (may be background-subtracted and scaled,
302 depending on configuration).
303 idFactory : `lsst.afw.table.IdFactory`
304 IdFactory to set source identifiers.
305 expId : `int`
306 Exposure identifier (integer) for RNG seed.
307 patchInfo : `lsst.skymap.PatchInfo`, optional
308 Description of the patch geometry. Only needed if
309 `~DetectCoaddSourceConfig.forceExactBinning` is `True`.
311 Returns
312 -------
313 result : `lsst.pipe.base.Struct`
314 Results as a struct with attributes:
316 ``sources``
317 Catalog of detections (`lsst.afw.table.SourceCatalog`).
318 ``backgrounds``
319 List of backgrounds (`list`).
320 """
321 if self.config.forceExactBinning:
322 exposure = self._cropToExactBinning(exposure, patchInfo)
323 if self.config.doScaleVariance:
324 varScale = self.scaleVariance.run(exposure.maskedImage)
325 exposure.getMetadata().add("VARIANCE_SCALE", varScale)
326 backgrounds = afwMath.BackgroundList()
327 table = afwTable.SourceTable.make(self.schema, idFactory)
328 detections = self.detection.run(table, exposure, expId=expId)
329 sources = detections.sources
330 if hasattr(detections, "background") and detections.background:
331 for bg in detections.background:
332 backgrounds.append(bg)
333 if len(backgrounds) == 0:
334 # Persist a constant background with value of NaN to get around
335 # inability to persist empty BackgroundList.
336 emptyBg = self._makeEmptyBackground(exposure, patchInfo)
337 backgrounds.append(emptyBg)
339 return Struct(outputSources=sources, outputBackgrounds=backgrounds, outputExposure=exposure)
341 def _cropToExactBinning(self, exposure, patchInfo):
342 """Crop a coadd `~lsst.afw.image.Exposure` instance to ensure exact
343 background binning.
345 Parameters
346 ----------
347 exposure : `lsst.afw.image.Exposure`
348 Exposure to crop, assumed to cover the patch outer bounding box.
349 patchInfo : `lsst.skymap.PatchInfo`
350 Description of the patch geometry.
352 Returns
353 -------
354 cropped : `lsst.afw.image.Exposure`
355 View of ``exposure`` with background bins that evenly divide both
356 the full cropped image and the patch inner region. The bounding
357 box is guaranteed to contain the patch inner bounding box and be
358 contained by the patch outer bounding box.
360 Raises
361 ------
362 ValueError
363 Raised if the patch inner region width or height is not a multiple
364 of the background bin size.
365 """
366 bbox = patchInfo.getInnerBBox()
367 if bbox.width % self.detection.background.binSizeX:
368 raise ValueError(
369 f"Patch inner width {bbox.width} does not evenly "
370 f"divide bin width {self.detection.background.binSizeX}."
371 )
372 if bbox.height % self.detection.background.binSizeY:
373 raise ValueError(
374 f"Patch inner height {bbox.height} does not evenly "
375 f"divide bin height {self.detection.background.binSizeY}."
376 )
377 outer_bbox = patchInfo.getOuterBBox()
378 n_bins_grow_x = (bbox.x.begin - outer_bbox.x.begin) // self.detection.background.binSizeX
379 n_bins_grow_y = (bbox.y.begin - outer_bbox.y.begin) // self.detection.background.binSizeY
380 bbox.grow(
381 Extent2I(
382 n_bins_grow_x*self.detection.background.binSizeX,
383 n_bins_grow_y*self.detection.background.binSizeY,
384 )
385 )
386 assert outer_bbox.contains(bbox)
387 assert bbox.contains(patchInfo.getInnerBBox())
388 assert bbox.width % self.detection.background.binSizeX == 0
389 assert bbox.height % self.detection.background.binSizeY == 0
390 return exposure[bbox]
392 def _makeEmptyBackground(self, exposure, patchInfo=None):
393 """Construct an empty `lsst.afw.math.BackgroundList` with NaN values.
395 Parameters
396 ----------
397 exposure : `lsst.afw.image.Exposure`
398 Exposure that the background should correspond to.
399 patchInfo : `lsst.skymap.PatchInfo`, optional
400 Description of the patch geometry. Only needed if
401 `~DetectCoaddSourceConfig.forceExactBinning` is `True`.
403 Returns
404 -------
405 background : `lsst.afw.math.BackgroundList`
406 A background object with a single layer and the same bin geometry
407 that a background for that exposure would have had if it had enough
408 usable pixels. This object cannot actually be used for background
409 subtraction.
410 """
411 # Create a backgroundList with one entry whose "stats image" is NaNs
412 # and has all pixels set as NO_DATA.
413 if self.config.forceExactBinning:
414 exposure = self._cropToExactBinning(exposure, patchInfo).clone()
416 bgLevel = np.nan
417 bgStats = afwImage.MaskedImageF(1, 1)
418 bgStats.set(bgLevel, 0, bgLevel)
419 bg = afwMath.BackgroundMI(exposure.getBBox(), bgStats)
420 bgData = (bg, afwMath.Interpolate.LINEAR, afwMath.REDUCE_INTERP_ORDER,
421 afwMath.ApproximateControl.UNKNOWN, 0, 0, False)
422 background = afwMath.BackgroundList()
423 background.append(bgData)
424 for bg, *_ in background:
425 stats = bg.getStatsImage()
426 stats.mask.array[:, :] = stats.mask.getPlaneBitMask("NO_DATA")
427 stats.variance.array[:, :] = 0.0
428 return background
431class MeasureMergedCoaddSourcesConnections(
432 PipelineTaskConnections,
433 dimensions=("tract", "patch", "band", "skymap"),
434 defaultTemplates={
435 "inputCoaddName": "deep",
436 "outputCoaddName": "deep",
437 "deblendedCatalog": "deblendedFlux",
438 },
439 deprecatedTemplates={
440 # TODO[DM-47797]: remove this deprecated connection template.
441 "deblendedCatalog": "Support for old deblender outputs will be removed after v29."
442 },
443):
444 inputSchema = cT.InitInput(
445 doc="Input schema for measure merged task produced by a deblender or detection task",
446 name="{inputCoaddName}Coadd_deblendedFlux_schema",
447 storageClass="SourceCatalog"
448 )
449 outputSchema = cT.InitOutput(
450 doc="Output schema after all new fields are added by task",
451 name="{inputCoaddName}Coadd_meas_schema",
452 storageClass="SourceCatalog"
453 )
454 # TODO[DM-47797]: remove this deprecated connection.
455 refCat = cT.PrerequisiteInput(
456 doc="Reference catalog used to match measured sources against known sources",
457 name="ref_cat",
458 storageClass="SimpleCatalog",
459 dimensions=("skypix",),
460 deferLoad=True,
461 multiple=True,
462 deprecated="Reference matching in measureCoaddSources will be removed after v29.",
463 )
464 exposure = cT.Input(
465 doc="Input non-cell-based coadd image",
466 name="{inputCoaddName}Coadd_calexp",
467 storageClass="ExposureF",
468 dimensions=("tract", "patch", "band", "skymap")
469 )
470 exposure_cells = cT.Input(
471 doc="Input cell-based coadd image",
472 name="{inputCoaddName}CoaddCell",
473 storageClass="MultipleCellCoadd",
474 dimensions=("tract", "patch", "band", "skymap"),
475 )
476 background = cT.Input(
477 doc="Background to subtract from cell-based coadd image",
478 name="{inputCoaddName}Coadd_calexp_background",
479 storageClass="Background",
480 dimensions=("tract", "patch", "band", "skymap")
481 )
482 skyMap = cT.Input(
483 doc="SkyMap to use in processing",
484 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
485 storageClass="SkyMap",
486 dimensions=("skymap",),
487 )
488 # TODO[DM-47424]: remove this deprecated connection.
489 visitCatalogs = cT.Input(
490 doc="Deprecated and unused.",
491 name="src",
492 dimensions=("instrument", "visit", "detector"),
493 storageClass="SourceCatalog",
494 multiple=True,
495 deprecated="Deprecated and unused. Will be removed after v29.",
496 )
497 sourceTableHandles = cT.Input(
498 doc=("Source tables that are derived from the ``CalibrateTask`` sources. "
499 "These tables contain astrometry and photometry flags, and optionally "
500 "PSF flags."),
501 name="sourceTable_visit",
502 storageClass="ArrowAstropy",
503 dimensions=("instrument", "visit"),
504 multiple=True,
505 deferLoad=True,
506 )
507 finalizedSourceTableHandles = cT.Input(
508 doc=("Finalized source tables from ``FinalizeCalibrationTask``. These "
509 "tables contain PSF flags from the finalized PSF estimation."),
510 name="finalized_src_table",
511 storageClass="ArrowAstropy",
512 dimensions=("instrument", "visit"),
513 multiple=True,
514 deferLoad=True,
515 )
516 finalVisitSummaryHandles = cT.Input(
517 doc="Final visit summary table",
518 name="finalVisitSummary",
519 storageClass="ExposureCatalog",
520 dimensions=("instrument", "visit"),
521 multiple=True,
522 deferLoad=True,
523 )
524 # TODO[DM-47797]: remove this deprecated connection.
525 inputCatalog = cT.Input(
526 doc=("Name of the input catalog to use."
527 "If the single band deblender was used this should be 'deblendedFlux."
528 "If the multi-band deblender was used this should be 'deblendedModel, "
529 "or deblendedFlux if the multiband deblender was configured to output "
530 "deblended flux catalogs. If no deblending was performed this should "
531 "be 'mergeDet'"),
532 name="{inputCoaddName}Coadd_{deblendedCatalog}",
533 storageClass="SourceCatalog",
534 deprecated="Support for old deblender outputs will be removed after v29.",
535 dimensions=("tract", "patch", "band", "skymap"),
536 )
537 scarletCatalog = cT.Input(
538 doc="Catalogs produced by multiband deblending",
539 name="{inputCoaddName}Coadd_deblendedCatalog",
540 storageClass="SourceCatalog",
541 dimensions=("tract", "patch", "skymap"),
542 )
543 scarletModels = cT.Input(
544 doc="Multiband scarlet models produced by the deblender",
545 name="{inputCoaddName}Coadd_scarletModelData",
546 storageClass="LsstScarletModelData",
547 dimensions=("tract", "patch", "skymap"),
548 )
549 outputSources = cT.Output(
550 doc="Source catalog containing all the measurement information generated in this task",
551 name="{outputCoaddName}Coadd_meas",
552 dimensions=("tract", "patch", "band", "skymap"),
553 storageClass="SourceCatalog",
554 )
555 # TODO[DM-47797]: remove this deprecated connection.
556 matchResult = cT.Output(
557 doc="Match catalog produced by configured matcher, optional on doMatchSources",
558 name="{outputCoaddName}Coadd_measMatch",
559 dimensions=("tract", "patch", "band", "skymap"),
560 storageClass="Catalog",
561 deprecated="Reference matching in measureCoaddSources will be removed after v29.",
562 )
563 # TODO[DM-47797]: remove this deprecated connection.
564 denormMatches = cT.Output(
565 doc="Denormalized Match catalog produced by configured matcher, optional on "
566 "doWriteMatchesDenormalized",
567 name="{outputCoaddName}Coadd_measMatchFull",
568 dimensions=("tract", "patch", "band", "skymap"),
569 storageClass="Catalog",
570 deprecated="Reference matching in measureCoaddSources will be removed after v29.",
571 )
573 def __init__(self, *, config=None):
574 super().__init__(config=config)
575 del self.visitCatalogs
576 if not config.doPropagateFlags:
577 del self.sourceTableHandles
578 del self.finalizedSourceTableHandles
579 del self.finalVisitSummaryHandles
580 else:
581 # Check for types of flags required.
582 if not config.propagateFlags.source_flags:
583 del self.sourceTableHandles
584 if not config.propagateFlags.finalized_source_flags:
585 del self.finalizedSourceTableHandles
586 # TODO[DM-47797]: only the 'if' block contents here should survive.
587 if config.inputCatalog == "deblendedCatalog":
588 del self.inputCatalog
589 if not config.doAddFootprints:
590 del self.scarletModels
591 else:
592 del self.deblendedCatalog
593 del self.scarletModels
595 # TODO[DM-47797]: delete the conditionals below.
596 if not config.doMatchSources:
597 del self.refCat
598 del self.matchResult
600 if not config.doWriteMatchesDenormalized:
601 del self.denormMatches
603 if config.useCellCoadds:
604 del self.exposure
605 else:
606 del self.exposure_cells
607 del self.background
610class MeasureMergedCoaddSourcesConfig(PipelineTaskConfig,
611 pipelineConnections=MeasureMergedCoaddSourcesConnections):
612 """Configuration parameters for the MeasureMergedCoaddSourcesTask
613 """
614 inputCatalog = ChoiceField(
615 dtype=str,
616 default="deblendedCatalog",
617 allowed={
618 "deblendedCatalog": "Output catalog from ScarletDeblendTask",
619 "deblendedFlux": "Output catalog from SourceDeblendTask",
620 "mergeDet": "The merged detections before deblending."
621 },
622 doc="The name of the input catalog.",
623 # TODO[DM-47797]: remove this config option and anything using it.
624 deprecated="Support for old deblender outputs will be removed after v29.",
625 )
626 doAddFootprints = Field(dtype=bool,
627 default=True,
628 doc="Whether or not to add footprints to the input catalog from scarlet models. "
629 "This should be true whenever using the multi-band deblender, "
630 "otherwise this should be False.")
631 doConserveFlux = Field(dtype=bool, default=True,
632 doc="Whether to use the deblender models as templates to re-distribute the flux "
633 "from the 'exposure' (True), or to perform measurements on the deblender "
634 "model footprints.")
635 doStripFootprints = Field(dtype=bool, default=True,
636 doc="Whether to strip footprints from the output catalog before "
637 "saving to disk. "
638 "This is usually done when using scarlet models to save disk space.")
639 useCellCoadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?")
640 measurement = ConfigurableField(target=SingleFrameMeasurementTask, doc="Source measurement")
641 setPrimaryFlags = ConfigurableField(target=SetPrimaryFlagsTask, doc="Set flags for primary tract/patch")
642 doPropagateFlags = Field(
643 dtype=bool, default=True,
644 doc="Whether to match sources to CCD catalogs to propagate flags (to e.g. identify PSF stars)"
645 )
646 propagateFlags = ConfigurableField(target=PropagateSourceFlagsTask, doc="Propagate source flags to coadd")
647 doMatchSources = Field(
648 dtype=bool,
649 default=False,
650 doc="Match sources to reference catalog?",
651 deprecated="Reference matching in measureCoaddSources will be removed after v29.",
652 )
653 match = ConfigurableField(
654 target=DirectMatchTask,
655 doc="Matching to reference catalog",
656 deprecated="Reference matching in measureCoaddSources will be removed after v29.",
657 )
658 doWriteMatchesDenormalized = Field(
659 dtype=bool,
660 default=False,
661 doc=("Write reference matches in denormalized format? "
662 "This format uses more disk space, but is more convenient to read."),
663 deprecated="Reference matching in measureCoaddSources will be removed after v29.",
664 )
665 coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
666 psfCache = Field(dtype=int, default=100, doc="Size of psfCache")
667 checkUnitsParseStrict = Field(
668 doc="Strictness of Astropy unit compatibility check, can be 'raise', 'warn' or 'silent'",
669 dtype=str,
670 default="raise",
671 )
672 doApCorr = Field(
673 dtype=bool,
674 default=True,
675 doc="Apply aperture corrections"
676 )
677 applyApCorr = ConfigurableField(
678 target=ApplyApCorrTask,
679 doc="Subtask to apply aperture corrections"
680 )
681 doRunCatalogCalculation = Field(
682 dtype=bool,
683 default=True,
684 doc='Run catalogCalculation task'
685 )
686 catalogCalculation = ConfigurableField(
687 target=CatalogCalculationTask,
688 doc="Subtask to run catalogCalculation plugins on catalog"
689 )
691 hasFakes = Field(
692 dtype=bool,
693 default=False,
694 doc="Should be set to True if fake sources have been inserted into the input data."
695 )
696 idGenerator = SkyMapIdGeneratorConfig.make_field()
698 @property
699 def refObjLoader(self):
700 return self.match.refObjLoader
702 def setDefaults(self):
703 super().setDefaults()
704 self.measurement.plugins.names |= ['base_InputCount',
705 'base_Variance',
706 'base_LocalPhotoCalib',
707 'base_LocalWcs']
709 # TODO: Remove STREAK in DM-44658, streak masking to happen only in
710 # ip_diffim; if we can propagate the streak mask from diffim, we can
711 # still set flags with it here.
712 self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['CLIPPED', 'SENSOR_EDGE',
713 'INEXACT_PSF']
714 self.measurement.plugins['base_PixelFlags'].masksFpCenter = ['CLIPPED', 'SENSOR_EDGE',
715 'INEXACT_PSF']
717 def validate(self):
718 super().validate()
720 if not self.doMatchSources and self.doWriteMatchesDenormalized:
721 raise ValueError("Cannot set doWriteMatchesDenormalized if doMatchSources is False.")
724class MeasureMergedCoaddSourcesTask(PipelineTask):
725 """Deblend sources from main catalog in each coadd seperately and measure.
727 Use peaks and footprints from a master catalog to perform deblending and
728 measurement in each coadd.
730 Given a master input catalog of sources (peaks and footprints) or deblender
731 outputs(including a HeavyFootprint in each band), measure each source on
732 the coadd. Repeating this procedure with the same master catalog across
733 multiple coadds will generate a consistent set of child sources.
735 The deblender retains all peaks and deblends any missing peaks (dropouts in
736 that band) as PSFs. Source properties are measured and the @c is-primary
737 flag (indicating sources with no children) is set. Visit flags are
738 propagated to the coadd sources.
740 Optionally, we can match the coadd sources to an external reference
741 catalog.
743 After MeasureMergedCoaddSourcesTask has been run on multiple coadds, we
744 have a set of per-band catalogs. The next stage in the multi-band
745 processing procedure will merge these measurements into a suitable catalog
746 for driving forced photometry.
748 Parameters
749 ----------
750 schema : ``lsst.afw.table.Schema`, optional
751 The schema of the merged detection catalog used as input to this one.
752 peakSchema : ``lsst.afw.table.Schema`, optional
753 The schema of the PeakRecords in the Footprints in the merged detection catalog.
754 refObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`, optional
755 An instance of ReferenceObjectLoader that supplies an external reference
756 catalog. May be None if the loader can be constructed from the butler argument or all steps
757 requiring a reference catalog are disabled.
758 initInputs : `dict`, optional
759 Dictionary that can contain a key ``inputSchema`` containing the
760 input schema. If present will override the value of ``schema``.
761 **kwargs
762 Additional keyword arguments.
763 """
765 _DefaultName = "measureCoaddSources"
766 ConfigClass = MeasureMergedCoaddSourcesConfig
768 def __init__(self, schema=None, peakSchema=None, refObjLoader=None, initInputs=None,
769 **kwargs):
770 super().__init__(**kwargs)
771 self.deblended = self.config.inputCatalog.startswith("deblended")
772 self.inputCatalog = "Coadd_" + self.config.inputCatalog
773 if initInputs is not None:
774 schema = initInputs['inputSchema'].schema
775 if schema is None:
776 raise ValueError("Schema must be defined.")
777 self.schemaMapper = afwTable.SchemaMapper(schema)
778 self.schemaMapper.addMinimalSchema(schema)
779 self.schema = self.schemaMapper.getOutputSchema()
780 self.algMetadata = PropertyList()
781 self.makeSubtask("measurement", schema=self.schema, algMetadata=self.algMetadata)
782 self.makeSubtask("setPrimaryFlags", schema=self.schema)
783 # TODO[DM-47797]: remove match subtask
784 if self.config.doMatchSources:
785 self.makeSubtask("match", refObjLoader=refObjLoader)
786 if self.config.doPropagateFlags:
787 self.makeSubtask("propagateFlags", schema=self.schema)
788 self.schema.checkUnits(parse_strict=self.config.checkUnitsParseStrict)
789 if self.config.doApCorr:
790 self.makeSubtask("applyApCorr", schema=self.schema)
791 if self.config.doRunCatalogCalculation:
792 self.makeSubtask("catalogCalculation", schema=self.schema)
794 self.outputSchema = afwTable.SourceCatalog(self.schema)
796 def runQuantum(self, butlerQC, inputRefs, outputRefs):
797 inputs = butlerQC.get(inputRefs)
799 # TODO[DM-47797]: remove this block
800 if self.config.doMatchSources:
801 refObjLoader = ReferenceObjectLoader([ref.datasetRef.dataId for ref in inputRefs.refCat],
802 inputs.pop('refCat'),
803 name=self.config.connections.refCat,
804 config=self.config.refObjLoader,
805 log=self.log)
806 self.match.setRefObjLoader(refObjLoader)
808 if self.config.useCellCoadds:
809 multiple_cell_coadd = inputs.pop("exposure_cells")
810 stitched_coadd = multiple_cell_coadd.stitch()
811 exposure = stitched_coadd.asExposure()
812 background = inputs.pop("background")
813 exposure.image -= background.getImage()
815 ccdInputs = stitched_coadd.ccds
816 apCorrMap = stitched_coadd.ap_corr_map
817 band = inputRefs.exposure_cells.dataId["band"]
818 else:
819 exposure = inputs.pop("exposure")
820 # Set psfcache
821 # move this to run after gen2 deprecation
822 exposure.getPsf().setCacheCapacity(self.config.psfCache)
824 ccdInputs = exposure.getInfo().getCoaddInputs().ccds
825 apCorrMap = exposure.getInfo().getApCorrMap()
826 band = inputRefs.exposure.dataId["band"]
828 # Get unique integer ID for IdFactory and RNG seeds; only the latter
829 # should really be used as the IDs all come from the input catalog.
830 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
832 # Transform inputCatalog
833 table = afwTable.SourceTable.make(self.schema, idGenerator.make_table_id_factory())
834 sources = afwTable.SourceCatalog(table)
835 # Load the correct input catalog
836 if "scarletCatalog" in inputs:
837 inputCatalog = inputs.pop("scarletCatalog")
838 catalogRef = inputRefs.scarletCatalog
839 else:
840 inputCatalog = inputs.pop("inputCatalog")
841 catalogRef = inputRefs.inputCatalog
842 sources.extend(inputCatalog, self.schemaMapper)
843 del inputCatalog
844 # Add the HeavyFootprints to the deblended sources
845 if self.config.doAddFootprints:
846 modelData = inputs.pop('scarletModels')
847 if self.config.doConserveFlux:
848 imageForRedistribution = exposure
849 else:
850 imageForRedistribution = None
851 updateCatalogFootprints(
852 modelData=modelData,
853 catalog=sources,
854 band=band,
855 imageForRedistribution=imageForRedistribution,
856 removeScarletData=True,
857 updateFluxColumns=True,
858 )
859 table = sources.getTable()
860 table.setMetadata(self.algMetadata) # Capture algorithm metadata to write out to the source catalog.
862 skyMap = inputs.pop('skyMap')
863 tractNumber = catalogRef.dataId['tract']
864 tractInfo = skyMap[tractNumber]
865 patchInfo = tractInfo.getPatchInfo(catalogRef.dataId['patch'])
866 skyInfo = Struct(
867 skyMap=skyMap,
868 tractInfo=tractInfo,
869 patchInfo=patchInfo,
870 wcs=tractInfo.getWcs(),
871 bbox=patchInfo.getOuterBBox()
872 )
874 sourceTableHandleDict = None
875 finalizedSourceTableHandleDict = None
876 finalVisitSummaryHandleDict = None
877 if self.config.doPropagateFlags:
878 if "sourceTableHandles" in inputs:
879 sourceTableHandles = inputs.pop("sourceTableHandles")
880 sourceTableHandleDict = {handle.dataId["visit"]: handle for handle in sourceTableHandles}
881 if "finalizedSourceTableHandles" in inputs:
882 finalizedSourceTableHandles = inputs.pop("finalizedSourceTableHandles")
883 finalizedSourceTableHandleDict = {handle.dataId["visit"]: handle
884 for handle in finalizedSourceTableHandles}
885 if "finalVisitSummaryHandles" in inputs:
886 finalVisitSummaryHandles = inputs.pop("finalVisitSummaryHandles")
887 finalVisitSummaryHandleDict = {handle.dataId["visit"]: handle
888 for handle in finalVisitSummaryHandles}
890 assert not inputs, "runQuantum got more inputs than expected."
891 outputs = self.run(
892 exposure=exposure,
893 sources=sources,
894 skyInfo=skyInfo,
895 exposureId=idGenerator.catalog_id,
896 ccdInputs=ccdInputs,
897 sourceTableHandleDict=sourceTableHandleDict,
898 finalizedSourceTableHandleDict=finalizedSourceTableHandleDict,
899 finalVisitSummaryHandleDict=finalVisitSummaryHandleDict,
900 apCorrMap=apCorrMap,
901 )
902 # Strip HeavyFootprints to save space on disk
903 if self.config.doStripFootprints:
904 sources = outputs.outputSources
905 for source in sources[sources["parent"] != 0]:
906 source.setFootprint(None)
907 butlerQC.put(outputs, outputRefs)
909 def run(self, exposure, sources, skyInfo, exposureId, ccdInputs=None,
910 sourceTableHandleDict=None, finalizedSourceTableHandleDict=None, finalVisitSummaryHandleDict=None,
911 apCorrMap=None):
912 """Run measurement algorithms on the input exposure, and optionally populate the
913 resulting catalog with extra information.
915 Parameters
916 ----------
917 exposure : `lsst.afw.exposure.Exposure`
918 The input exposure on which measurements are to be performed.
919 sources : `lsst.afw.table.SourceCatalog`
920 A catalog built from the results of merged detections, or
921 deblender outputs.
922 parentCatalog : `lsst.afw.table.SourceCatalog`
923 Catalog of parent sources corresponding to sources.
924 skyInfo : `lsst.pipe.base.Struct`
925 A struct containing information about the position of the input exposure within
926 a `SkyMap`, the `SkyMap`, its `Wcs`, and its bounding box.
927 exposureId : `int` or `bytes`
928 Packed unique number or bytes unique to the input exposure.
929 ccdInputs : `lsst.afw.table.ExposureCatalog`, optional
930 Catalog containing information on the individual visits which went into making
931 the coadd.
932 sourceTableHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
933 Dict for sourceTable_visit handles (key is visit) for propagating flags.
934 These tables contain astrometry and photometry flags, and optionally PSF flags.
935 finalizedSourceTableHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
936 Dict for finalized_src_table handles (key is visit) for propagating flags.
937 These tables contain PSF flags from the finalized PSF estimation.
938 finalVisitSummaryHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
939 Dict for visit_summary handles (key is visit) for visit-level information.
940 These tables contain the WCS information of the single-visit input images.
941 apCorrMap : `lsst.afw.image.ApCorrMap`, optional
942 Aperture correction map attached to the ``exposure``. If None, it
943 will be read from the ``exposure``.
945 Returns
946 -------
947 results : `lsst.pipe.base.Struct`
948 Results of running measurement task. Will contain the catalog in the
949 sources attribute. Optionally will have results of matching to a
950 reference catalog in the matchResults attribute, and denormalized
951 matches in the denormMatches attribute.
952 """
953 if self.config.doPropagateFlags:
954 # These mask planes may not be defined on the coadds always.
955 # We add the mask planes, which is a no-op if already defined.
956 for maskPlane in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
957 exposure.mask.addMaskPlane(maskPlane)
958 for maskPlane in self.config.measurement.plugins["base_PixelFlags"].masksFpCenter:
959 exposure.mask.addMaskPlane(maskPlane)
961 self.measurement.run(sources, exposure, exposureId=exposureId)
963 if self.config.doApCorr:
964 if apCorrMap is None:
965 apCorrMap = exposure.getInfo().getApCorrMap()
966 self.applyApCorr.run(
967 catalog=sources,
968 apCorrMap=apCorrMap,
969 )
971 # TODO DM-11568: this contiguous check-and-copy could go away if we
972 # reserve enough space during SourceDetection and/or SourceDeblend.
973 # NOTE: sourceSelectors require contiguous catalogs, so ensure
974 # contiguity now, so views are preserved from here on.
975 if not sources.isContiguous():
976 sources = sources.copy(deep=True)
978 if self.config.doRunCatalogCalculation:
979 self.catalogCalculation.run(sources)
981 self.setPrimaryFlags.run(sources, skyMap=skyInfo.skyMap, tractInfo=skyInfo.tractInfo,
982 patchInfo=skyInfo.patchInfo)
983 if self.config.doPropagateFlags:
984 self.propagateFlags.run(
985 sources,
986 ccdInputs,
987 sourceTableHandleDict,
988 finalizedSourceTableHandleDict,
989 finalVisitSummaryHandleDict,
990 )
992 results = Struct()
994 # TODO[DM-47797]: remove this block
995 if self.config.doMatchSources:
996 matchResult = self.match.run(sources, exposure.getInfo().getFilter().bandLabel)
997 matches = afwTable.packMatches(matchResult.matches)
998 matches.table.setMetadata(matchResult.matchMeta)
999 results.matchResult = matches
1000 if self.config.doWriteMatchesDenormalized:
1001 if matchResult.matches:
1002 denormMatches = denormalizeMatches(matchResult.matches, matchResult.matchMeta)
1003 else:
1004 self.log.warning("No matches, so generating dummy denormalized matches file")
1005 denormMatches = afwTable.BaseCatalog(afwTable.Schema())
1006 denormMatches.setMetadata(PropertyList())
1007 denormMatches.getMetadata().add("COMMENT",
1008 "This catalog is empty because no matches were found.")
1009 results.denormMatches = denormMatches
1010 results.denormMatches = denormMatches
1012 results.outputSources = sources
1013 return results