Coverage for python/lsst/pipe/tasks/multiBand.py: 23%
298 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-25 09:10 +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
38from lsst.meas.algorithms import (
39 DynamicDetectionTask,
40 ExceedsMaxVarianceScaleError,
41 InsufficientSourcesError,
42 PsfGenerationError,
43 ScaleVarianceTask,
44 SetPrimaryFlagsTask,
45 TooManyMaskedPixelsError,
46 ZeroFootprintError,
47)
48from lsst.meas.base import (
49 SingleFrameMeasurementTask,
50 ApplyApCorrTask,
51 CatalogCalculationTask,
52 SkyMapIdGeneratorConfig,
53)
54from lsst.meas.extensions.scarlet.io import updateCatalogFootprints
55from lsst.pipe.tasks.propagateSourceFlags import PropagateSourceFlagsTask
56import lsst.afw.image as afwImage
57import lsst.afw.math as afwMath
58import lsst.afw.table as afwTable
59from lsst.daf.base import PropertyList
60from lsst.skymap import BaseSkyMap
62# NOTE: these imports are a convenience so multiband users only have to import this file.
63from .mergeDetections import MergeDetectionsConfig, MergeDetectionsTask # noqa: F401
64from .mergeMeasurements import MergeMeasurementsConfig, MergeMeasurementsTask # noqa: F401
65from .multiBandUtils import CullPeaksConfig # noqa: F401
66from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiConfig # noqa: F401
67from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiTask # noqa: F401
70"""
71New set types:
72* deepCoadd_det: detections from what used to be processCoadd (tract, patch, filter)
73* deepCoadd_mergeDet: merged detections (tract, patch)
74* deepCoadd_meas: measurements of merged detections (tract, patch, filter)
75* deepCoadd_ref: reference sources (tract, patch)
76All of these have associated *_schema catalogs that require no data ID and hold no records.
78In addition, we have a schema-only dataset, which saves the schema for the PeakRecords in
79the mergeDet, meas, and ref dataset Footprints:
80* deepCoadd_peak_schema
81"""
84##############################################################################################################
85class DetectCoaddSourcesConnections(PipelineTaskConnections,
86 dimensions=("tract", "patch", "band", "skymap"),
87 defaultTemplates={"inputCoaddName": "deep", "outputCoaddName": "deep"}):
88 detectionSchema = cT.InitOutput(
89 doc="Schema of the detection catalog",
90 name="{outputCoaddName}Coadd_det_schema",
91 storageClass="SourceCatalog",
92 )
93 exposure = cT.Input(
94 doc="Exposure on which detections are to be performed. ",
95 name="{inputCoaddName}Coadd",
96 storageClass="ExposureF",
97 dimensions=("tract", "patch", "band", "skymap")
98 )
99 exposure_cells = cT.Input(
100 doc="Exposure on which detections are to be performed. ",
101 name="{inputCoaddName}CoaddCell",
102 storageClass="MultipleCellCoadd",
103 dimensions=("tract", "patch", "band", "skymap"),
104 )
105 skyMap = cT.Input(
106 doc="Description of the skymap's tracts and patches.",
107 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
108 storageClass="SkyMap",
109 dimensions=("skymap",),
110 )
111 outputBackgrounds = cT.Output(
112 doc="Output Backgrounds used in detection",
113 name="{outputCoaddName}Coadd_calexp_background",
114 storageClass="Background",
115 dimensions=("tract", "patch", "band", "skymap")
116 )
117 outputSources = cT.Output(
118 doc="Detected sources catalog",
119 name="{outputCoaddName}Coadd_det",
120 storageClass="SourceCatalog",
121 dimensions=("tract", "patch", "band", "skymap")
122 )
123 outputExposure = cT.Output(
124 doc="Exposure post detection",
125 name="{outputCoaddName}Coadd_calexp",
126 storageClass="ExposureF",
127 dimensions=("tract", "patch", "band", "skymap")
128 )
130 def __init__(self, *, config=None):
131 super().__init__(config=config)
132 assert isinstance(config, DetectCoaddSourcesConfig)
134 if config.useCellCoadds:
135 del self.exposure
136 else:
137 del self.exposure_cells
139 if not self.config.forceExactBinning:
140 del self.skyMap
141 if self.config.writeOnlyBackgrounds:
142 del self.outputExposure
143 del self.outputSources
144 del self.detectionSchema
147class DetectCoaddSourcesConfig(PipelineTaskConfig, pipelineConnections=DetectCoaddSourcesConnections):
148 """Configuration parameters for the DetectCoaddSourcesTask
149 """
151 doScaleVariance = Field(dtype=bool, default=True, doc="Scale variance plane using empirical noise?")
152 scaleVariance = ConfigurableField(target=ScaleVarianceTask, doc="Variance rescaling")
153 detection = ConfigurableField(target=DynamicDetectionTask, doc="Source detection")
154 coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
155 useCellCoadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?")
156 hasFakes = Field(
157 dtype=bool,
158 default=False,
159 doc="Should be set to True if fake sources have been inserted into the input data.",
160 )
161 idGenerator = SkyMapIdGeneratorConfig.make_field()
162 forceExactBinning = Field(
163 dtype=bool,
164 default=False,
165 doc=(
166 "Check that the background bin size evenly divides the patch inner region, and "
167 "crop the outer region to an integer number of bins."
168 )
169 )
170 writeOnlyBackgrounds = Field(dtype=bool, default=False, doc="If true, only save the background models.")
171 writeEmptyBackgrounds = Field(
172 dtype=bool,
173 default=True,
174 doc=(
175 "If true, save a placeholder background with NaNs in all bins (but the right geometry) when "
176 "there are no pixels to compute a background from. This can be useful if a later task combines "
177 "backgrounds from multiple patches as input."
178 )
179 )
181 def setDefaults(self):
182 super().setDefaults()
183 self.detection.thresholdType = "pixel_stdev"
184 self.detection.isotropicGrow = True
185 # Coadds are made from background-subtracted CCDs, so any background subtraction should be very basic
186 self.detection.reEstimateBackground = False
187 self.detection.background.useApprox = False
188 self.detection.background.binSize = 4096
189 self.detection.background.undersampleStyle = 'REDUCE_INTERP_ORDER'
190 self.detection.doTempWideBackground = True # Suppress large footprints that overwhelm the deblender
191 # Include band in packed data IDs that go into object IDs (None -> "as
192 # many bands as are defined", rather than the default of zero).
193 self.idGenerator.packer.n_bands = None
196class DetectCoaddSourcesTask(PipelineTask):
197 """Detect sources on a single filter coadd.
199 Coadding individual visits requires each exposure to be warped. This
200 introduces covariance in the noise properties across pixels. Before
201 detection, we correct the coadd variance by scaling the variance plane in
202 the coadd to match the observed variance. This is an approximate
203 approach -- strictly, we should propagate the full covariance matrix --
204 but it is simple and works well in practice.
206 After scaling the variance plane, we detect sources and generate footprints
207 by delegating to the @ref SourceDetectionTask_ "detection" subtask.
209 DetectCoaddSourcesTask is meant to be run after assembling a coadded image
210 in a given band. The purpose of the task is to update the background,
211 detect all sources in a single band and generate a set of parent
212 footprints. Subsequent tasks in the multi-band processing procedure will
213 merge sources across bands and, eventually, perform forced photometry.
215 Parameters
216 ----------
217 schema : `lsst.afw.table.Schema`, optional
218 Initial schema for the output catalog, modified-in place to include all
219 fields set by this task. If None, the source minimal schema will be used.
220 **kwargs
221 Additional keyword arguments.
222 """
224 _DefaultName = "detectCoaddSources"
225 ConfigClass = DetectCoaddSourcesConfig
227 def __init__(self, schema=None, **kwargs):
228 # N.B. Super is used here to handle the multiple inheritance of PipelineTasks, the init tree
229 # call structure has been reviewed carefully to be sure super will work as intended.
230 super().__init__(**kwargs)
231 if schema is None:
232 schema = afwTable.SourceTable.makeMinimalSchema()
233 self.schema = schema
234 self.makeSubtask("detection", schema=self.schema)
235 if self.config.doScaleVariance:
236 self.makeSubtask("scaleVariance")
238 self.detectionSchema = afwTable.SourceCatalog(self.schema)
240 def runQuantum(self, butlerQC, inputRefs, outputRefs):
241 inputs = butlerQC.get(inputRefs)
242 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
244 if self.config.useCellCoadds:
245 multiple_cell_coadd = inputs.pop("exposure_cells")
246 exposure = multiple_cell_coadd.stitch().asExposure()
247 else:
248 exposure = inputs.pop("exposure")
250 skyMap = inputs.pop("skyMap", None)
251 if skyMap is not None:
252 patchInfo = skyMap[butlerQC.quantum.dataId["tract"]][butlerQC.quantum.dataId["patch"]]
253 else:
254 patchInfo = None
256 assert not inputs, "runQuantum got more inputs than expected."
257 try:
258 outputs = self.run(
259 exposure=exposure,
260 idFactory=idGenerator.make_table_id_factory(),
261 expId=idGenerator.catalog_id,
262 patchInfo=patchInfo,
263 )
264 except (
265 TooManyMaskedPixelsError,
266 ExceedsMaxVarianceScaleError,
267 InsufficientSourcesError,
268 PsfGenerationError,
269 ZeroFootprintError,
270 ) as e:
271 if self.config.writeEmptyBackgrounds:
272 butlerQC.put(self._makeEmptyBackground(exposure, patchInfo), outputRefs.outputBackgrounds)
273 # Detection failed, so clear any leftover the detected mask planes.
274 for maskName in ["DETECTED", "DETECTED_NEGATIVE"]:
275 if maskName in exposure.mask.getMaskPlaneDict().keys():
276 detectedMask = exposure.mask.getMaskPlane(maskName)
277 exposure.mask.clearMaskPlane(detectedMask)
278 butlerQC.put(exposure, outputRefs.outputExposure)
279 error = AnnotatedPartialOutputsError.annotate(
280 e,
281 self,
282 exposure,
283 log=self.log,
284 )
285 raise error from e
287 butlerQC.put(outputs, outputRefs)
289 def run(self, exposure, idFactory, expId, patchInfo=None):
290 """Run detection on an exposure.
292 First scale the variance plane to match the observed variance
293 using ``ScaleVarianceTask``. Then invoke the ``SourceDetectionTask_`` "detection" subtask to
294 detect sources.
296 Parameters
297 ----------
298 exposure : `lsst.afw.image.Exposure`
299 Exposure on which to detect (may be background-subtracted and scaled,
300 depending on configuration).
301 idFactory : `lsst.afw.table.IdFactory`
302 IdFactory to set source identifiers.
303 expId : `int`
304 Exposure identifier (integer) for RNG seed.
305 patchInfo : `lsst.skymap.PatchInfo`, optional
306 Description of the patch geometry. Only needed if
307 `~DetectCoaddSourceConfig.forceExactBinning` is `True`.
309 Returns
310 -------
311 result : `lsst.pipe.base.Struct`
312 Results as a struct with attributes:
314 ``sources``
315 Catalog of detections (`lsst.afw.table.SourceCatalog`).
316 ``backgrounds``
317 List of backgrounds (`list`).
318 """
319 if self.config.forceExactBinning:
320 exposure = self._cropToExactBinning(exposure, patchInfo)
321 if self.config.doScaleVariance:
322 varScale = self.scaleVariance.run(exposure.maskedImage)
323 exposure.getMetadata().add("VARIANCE_SCALE", varScale)
324 backgrounds = afwMath.BackgroundList()
325 table = afwTable.SourceTable.make(self.schema, idFactory)
326 detections = self.detection.run(table, exposure, expId=expId)
327 sources = detections.sources
328 if hasattr(detections, "background") and detections.background:
329 for bg in detections.background:
330 backgrounds.append(bg)
331 if len(backgrounds) == 0:
332 # Persist a constant background with value of NaN to get around
333 # inability to persist empty BackgroundList.
334 emptyBg = self._makeEmptyBackground(exposure, patchInfo)
335 backgrounds.append(emptyBg)
337 return Struct(outputSources=sources, outputBackgrounds=backgrounds, outputExposure=exposure)
339 def _cropToExactBinning(self, exposure, patchInfo):
340 """Crop a coadd `~lsst.afw.image.Exposure` instance to ensure exact
341 background binning.
343 Parameters
344 ----------
345 exposure : `lsst.afw.image.Exposure`
346 Exposure to crop, assumed to cover the patch outer bounding box.
347 patchInfo : `lsst.skymap.PatchInfo`
348 Description of the patch geometry.
350 Returns
351 -------
352 cropped : `lsst.afw.image.Exposure`
353 View of ``exposure`` with background bins that evenly divide both
354 the full cropped image and the patch inner region. The bounding
355 box is guaranteed to contain the patch inner bounding box and be
356 contained by the patch outer bounding box.
358 Raises
359 ------
360 ValueError
361 Raised if the patch inner region width or height is not a multiple
362 of the background bin size.
363 """
364 bbox = patchInfo.getInnerBBox()
365 if bbox.width % self.detection.background.binSizeX:
366 raise ValueError(
367 f"Patch inner width {bbox.width} does not evenly "
368 f"divide bin width {self.detection.background.binSizeX}."
369 )
370 if bbox.height % self.detection.background.binSizeY:
371 raise ValueError(
372 f"Patch inner height {bbox.height} does not evenly "
373 f"divide bin height {self.detection.background.binSizeY}."
374 )
375 outer_bbox = patchInfo.getOuterBBox()
376 n_bins_grow_x = (bbox.x.begin - outer_bbox.x.begin) // self.detection.background.binSizeX
377 n_bins_grow_y = (bbox.y.begin - outer_bbox.y.begin) // self.detection.background.binSizeY
378 bbox.grow(
379 Extent2I(
380 n_bins_grow_x*self.detection.background.binSizeX,
381 n_bins_grow_y*self.detection.background.binSizeY,
382 )
383 )
384 assert outer_bbox.contains(bbox)
385 assert bbox.contains(patchInfo.getInnerBBox())
386 assert bbox.width % self.detection.background.binSizeX == 0
387 assert bbox.height % self.detection.background.binSizeY == 0
388 return exposure[bbox]
390 def _makeEmptyBackground(self, exposure, patchInfo=None):
391 """Construct an empty `lsst.afw.math.BackgroundList` with NaN values.
393 Parameters
394 ----------
395 exposure : `lsst.afw.image.Exposure`
396 Exposure that the background should correspond to.
397 patchInfo : `lsst.skymap.PatchInfo`, optional
398 Description of the patch geometry. Only needed if
399 `~DetectCoaddSourceConfig.forceExactBinning` is `True`.
401 Returns
402 -------
403 background : `lsst.afw.math.BackgroundList`
404 A background object with a single layer and the same bin geometry
405 that a background for that exposure would have had if it had enough
406 usable pixels. This object cannot actually be used for background
407 subtraction.
408 """
409 # Create a backgroundList with one entry whose "stats image" is NaNs
410 # and has all pixels set as NO_DATA.
411 if self.config.forceExactBinning:
412 exposure = self._cropToExactBinning(exposure, patchInfo).clone()
414 bgLevel = np.nan
415 bgStats = afwImage.MaskedImageF(1, 1)
416 bgStats.set(bgLevel, 0, bgLevel)
417 bg = afwMath.BackgroundMI(exposure.getBBox(), bgStats)
418 bgData = (bg, afwMath.Interpolate.LINEAR, afwMath.REDUCE_INTERP_ORDER,
419 afwMath.ApproximateControl.UNKNOWN, 0, 0, False)
420 background = afwMath.BackgroundList()
421 background.append(bgData)
422 for bg, *_ in background:
423 stats = bg.getStatsImage()
424 stats.mask.array[:, :] = stats.mask.getPlaneBitMask("NO_DATA")
425 stats.variance.array[:, :] = 0.0
426 return background
429class MeasureMergedCoaddSourcesConnections(
430 PipelineTaskConnections,
431 dimensions=("tract", "patch", "band", "skymap"),
432 defaultTemplates={
433 "inputCoaddName": "deep",
434 "outputCoaddName": "deep",
435 },
436):
437 inputSchema = cT.InitInput(
438 doc="Input schema for measure merged task produced by a deblender or detection task",
439 name="{inputCoaddName}Coadd_deblendedFlux_schema",
440 storageClass="SourceCatalog"
441 )
442 outputSchema = cT.InitOutput(
443 doc="Output schema after all new fields are added by task",
444 name="{inputCoaddName}Coadd_meas_schema",
445 storageClass="SourceCatalog"
446 )
447 exposure = cT.Input(
448 doc="Input non-cell-based coadd image",
449 name="{inputCoaddName}Coadd_calexp",
450 storageClass="ExposureF",
451 dimensions=("tract", "patch", "band", "skymap")
452 )
453 exposure_cells = cT.Input(
454 doc="Input cell-based coadd image",
455 name="{inputCoaddName}CoaddCell",
456 storageClass="MultipleCellCoadd",
457 dimensions=("tract", "patch", "band", "skymap"),
458 )
459 background = cT.Input(
460 doc="Background to subtract from cell-based coadd image",
461 name="{inputCoaddName}Coadd_calexp_background",
462 storageClass="Background",
463 dimensions=("tract", "patch", "band", "skymap")
464 )
465 skyMap = cT.Input(
466 doc="SkyMap to use in processing",
467 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
468 storageClass="SkyMap",
469 dimensions=("skymap",),
470 )
471 sourceTableHandles = cT.Input(
472 doc=("Source tables that are derived from the ``CalibrateTask`` sources. "
473 "These tables contain astrometry and photometry flags, and optionally "
474 "PSF flags."),
475 name="sourceTable_visit",
476 storageClass="ArrowAstropy",
477 dimensions=("instrument", "visit"),
478 multiple=True,
479 deferLoad=True,
480 )
481 finalizedSourceTableHandles = cT.Input(
482 doc=("Finalized source tables from ``FinalizeCalibrationTask``. These "
483 "tables contain PSF flags from the finalized PSF estimation."),
484 name="finalized_src_table",
485 storageClass="ArrowAstropy",
486 dimensions=("instrument", "visit"),
487 multiple=True,
488 deferLoad=True,
489 )
490 finalVisitSummaryHandles = cT.Input(
491 doc="Final visit summary table",
492 name="finalVisitSummary",
493 storageClass="ExposureCatalog",
494 dimensions=("instrument", "visit"),
495 multiple=True,
496 deferLoad=True,
497 )
498 scarletCatalog = cT.Input(
499 doc="Catalogs produced by multiband deblending",
500 name="{inputCoaddName}Coadd_deblendedCatalog",
501 storageClass="SourceCatalog",
502 dimensions=("tract", "patch", "skymap"),
503 )
504 scarletModels = cT.Input(
505 doc="Multiband scarlet models produced by the deblender",
506 name="{inputCoaddName}Coadd_scarletModelData",
507 storageClass="LsstScarletModelData",
508 dimensions=("tract", "patch", "skymap"),
509 )
510 outputSources = cT.Output(
511 doc="Source catalog containing all the measurement information generated in this task",
512 name="{outputCoaddName}Coadd_meas",
513 dimensions=("tract", "patch", "band", "skymap"),
514 storageClass="SourceCatalog",
515 )
517 def __init__(self, *, config=None):
518 super().__init__(config=config)
519 if not config.doPropagateFlags:
520 del self.sourceTableHandles
521 del self.finalizedSourceTableHandles
522 del self.finalVisitSummaryHandles
523 else:
524 # Check for types of flags required.
525 if not config.propagateFlags.source_flags:
526 del self.sourceTableHandles
527 if not config.propagateFlags.finalized_source_flags:
528 del self.finalizedSourceTableHandles
529 if not config.doAddFootprints:
530 del self.scarletModels
532 if config.useCellCoadds:
533 del self.exposure
534 else:
535 del self.exposure_cells
536 del self.background
539class MeasureMergedCoaddSourcesConfig(PipelineTaskConfig,
540 pipelineConnections=MeasureMergedCoaddSourcesConnections):
541 """Configuration parameters for the MeasureMergedCoaddSourcesTask
542 """
543 doAddFootprints = Field(dtype=bool,
544 default=True,
545 doc="Whether or not to add footprints to the input catalog from scarlet models. "
546 "This should be true whenever using the multi-band deblender, "
547 "otherwise this should be False.")
548 doConserveFlux = Field(dtype=bool, default=True,
549 doc="Whether to use the deblender models as templates to re-distribute the flux "
550 "from the 'exposure' (True), or to perform measurements on the deblender "
551 "model footprints.")
552 doStripFootprints = Field(dtype=bool, default=True,
553 doc="Whether to strip footprints from the output catalog before "
554 "saving to disk. "
555 "This is usually done when using scarlet models to save disk space.")
556 useCellCoadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?")
557 measurement = ConfigurableField(target=SingleFrameMeasurementTask, doc="Source measurement")
558 setPrimaryFlags = ConfigurableField(target=SetPrimaryFlagsTask, doc="Set flags for primary tract/patch")
559 doPropagateFlags = Field(
560 dtype=bool, default=True,
561 doc="Whether to match sources to CCD catalogs to propagate flags (to e.g. identify PSF stars)"
562 )
563 propagateFlags = ConfigurableField(target=PropagateSourceFlagsTask, doc="Propagate source flags to coadd")
564 coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
565 psfCache = Field(dtype=int, default=100, doc="Size of psfCache")
566 checkUnitsParseStrict = Field(
567 doc="Strictness of Astropy unit compatibility check, can be 'raise', 'warn' or 'silent'",
568 dtype=str,
569 default="raise",
570 )
571 doApCorr = Field(
572 dtype=bool,
573 default=True,
574 doc="Apply aperture corrections"
575 )
576 applyApCorr = ConfigurableField(
577 target=ApplyApCorrTask,
578 doc="Subtask to apply aperture corrections"
579 )
580 doRunCatalogCalculation = Field(
581 dtype=bool,
582 default=True,
583 doc='Run catalogCalculation task'
584 )
585 catalogCalculation = ConfigurableField(
586 target=CatalogCalculationTask,
587 doc="Subtask to run catalogCalculation plugins on catalog"
588 )
590 hasFakes = Field(
591 dtype=bool,
592 default=False,
593 doc="Should be set to True if fake sources have been inserted into the input data."
594 )
595 idGenerator = SkyMapIdGeneratorConfig.make_field()
597 def setDefaults(self):
598 super().setDefaults()
599 self.measurement.plugins.names |= ['base_InputCount',
600 'base_Variance',
601 'base_LocalPhotoCalib',
602 'base_LocalWcs']
604 # TODO: Remove STREAK in DM-44658, streak masking to happen only in
605 # ip_diffim; if we can propagate the streak mask from diffim, we can
606 # still set flags with it here.
607 self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['CLIPPED', 'SENSOR_EDGE',
608 'INEXACT_PSF']
609 self.measurement.plugins['base_PixelFlags'].masksFpCenter = ['CLIPPED', 'SENSOR_EDGE',
610 'INEXACT_PSF']
613class MeasureMergedCoaddSourcesTask(PipelineTask):
614 """Deblend sources from main catalog in each coadd seperately and measure.
616 Use peaks and footprints from a master catalog to perform deblending and
617 measurement in each coadd.
619 Given a master input catalog of sources (peaks and footprints) or deblender
620 outputs(including a HeavyFootprint in each band), measure each source on
621 the coadd. Repeating this procedure with the same master catalog across
622 multiple coadds will generate a consistent set of child sources.
624 The deblender retains all peaks and deblends any missing peaks (dropouts in
625 that band) as PSFs. Source properties are measured and the @c is-primary
626 flag (indicating sources with no children) is set. Visit flags are
627 propagated to the coadd sources.
629 After MeasureMergedCoaddSourcesTask has been run on multiple coadds, we
630 have a set of per-band catalogs. The next stage in the multi-band
631 processing procedure will merge these measurements into a suitable catalog
632 for driving forced photometry.
634 Parameters
635 ----------
636 schema : ``lsst.afw.table.Schema`, optional
637 The schema of the merged detection catalog used as input to this one.
638 peakSchema : ``lsst.afw.table.Schema`, optional
639 The schema of the PeakRecords in the Footprints in the merged detection catalog.
640 initInputs : `dict`, optional
641 Dictionary that can contain a key ``inputSchema`` containing the
642 input schema. If present will override the value of ``schema``.
643 **kwargs
644 Additional keyword arguments.
645 """
647 _DefaultName = "measureCoaddSources"
648 ConfigClass = MeasureMergedCoaddSourcesConfig
650 def __init__(self, schema=None, peakSchema=None, initInputs=None, **kwargs):
651 super().__init__(**kwargs)
652 if initInputs is not None:
653 schema = initInputs['inputSchema'].schema
654 if schema is None:
655 raise ValueError("Schema must be defined.")
656 self.schemaMapper = afwTable.SchemaMapper(schema)
657 self.schemaMapper.addMinimalSchema(schema)
658 self.schema = self.schemaMapper.getOutputSchema()
659 self.algMetadata = PropertyList()
660 self.makeSubtask("measurement", schema=self.schema, algMetadata=self.algMetadata)
661 self.makeSubtask("setPrimaryFlags", schema=self.schema)
662 if self.config.doPropagateFlags:
663 self.makeSubtask("propagateFlags", schema=self.schema)
664 self.schema.checkUnits(parse_strict=self.config.checkUnitsParseStrict)
665 if self.config.doApCorr:
666 self.makeSubtask("applyApCorr", schema=self.schema)
667 if self.config.doRunCatalogCalculation:
668 self.makeSubtask("catalogCalculation", schema=self.schema)
670 self.outputSchema = afwTable.SourceCatalog(self.schema)
672 def runQuantum(self, butlerQC, inputRefs, outputRefs):
673 inputs = butlerQC.get(inputRefs)
675 if self.config.useCellCoadds:
676 multiple_cell_coadd = inputs.pop("exposure_cells")
677 stitched_coadd = multiple_cell_coadd.stitch()
678 exposure = stitched_coadd.asExposure()
679 background = inputs.pop("background")
680 exposure.image -= background.getImage()
682 ccdInputs = stitched_coadd.ccds
683 apCorrMap = stitched_coadd.ap_corr_map
684 band = inputRefs.exposure_cells.dataId["band"]
685 else:
686 exposure = inputs.pop("exposure")
687 # Set psfcache
688 # move this to run after gen2 deprecation
689 exposure.getPsf().setCacheCapacity(self.config.psfCache)
691 ccdInputs = exposure.getInfo().getCoaddInputs().ccds
692 apCorrMap = exposure.getInfo().getApCorrMap()
693 band = inputRefs.exposure.dataId["band"]
695 # Get unique integer ID for IdFactory and RNG seeds; only the latter
696 # should really be used as the IDs all come from the input catalog.
697 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
699 # Transform inputCatalog
700 table = afwTable.SourceTable.make(self.schema, idGenerator.make_table_id_factory())
701 sources = afwTable.SourceCatalog(table)
702 # Load the correct input catalog
703 inputCatalog = inputs.pop("scarletCatalog")
704 catalogRef = inputRefs.scarletCatalog
705 sources.extend(inputCatalog, self.schemaMapper)
706 del inputCatalog
707 # Add the HeavyFootprints to the deblended sources
708 if self.config.doAddFootprints:
709 modelData = inputs.pop('scarletModels')
710 if self.config.doConserveFlux:
711 imageForRedistribution = exposure
712 else:
713 imageForRedistribution = None
714 updateCatalogFootprints(
715 modelData=modelData,
716 catalog=sources,
717 band=band,
718 imageForRedistribution=imageForRedistribution,
719 removeScarletData=True,
720 updateFluxColumns=True,
721 )
722 table = sources.getTable()
723 table.setMetadata(self.algMetadata) # Capture algorithm metadata to write out to the source catalog.
725 skyMap = inputs.pop('skyMap')
726 tractNumber = catalogRef.dataId['tract']
727 tractInfo = skyMap[tractNumber]
728 patchInfo = tractInfo.getPatchInfo(catalogRef.dataId['patch'])
729 skyInfo = Struct(
730 skyMap=skyMap,
731 tractInfo=tractInfo,
732 patchInfo=patchInfo,
733 wcs=tractInfo.getWcs(),
734 bbox=patchInfo.getOuterBBox()
735 )
737 sourceTableHandleDict = None
738 finalizedSourceTableHandleDict = None
739 finalVisitSummaryHandleDict = None
740 if self.config.doPropagateFlags:
741 if "sourceTableHandles" in inputs:
742 sourceTableHandles = inputs.pop("sourceTableHandles")
743 sourceTableHandleDict = {handle.dataId["visit"]: handle for handle in sourceTableHandles}
744 if "finalizedSourceTableHandles" in inputs:
745 finalizedSourceTableHandles = inputs.pop("finalizedSourceTableHandles")
746 finalizedSourceTableHandleDict = {handle.dataId["visit"]: handle
747 for handle in finalizedSourceTableHandles}
748 if "finalVisitSummaryHandles" in inputs:
749 finalVisitSummaryHandles = inputs.pop("finalVisitSummaryHandles")
750 finalVisitSummaryHandleDict = {handle.dataId["visit"]: handle
751 for handle in finalVisitSummaryHandles}
753 assert not inputs, "runQuantum got more inputs than expected."
754 outputs = self.run(
755 exposure=exposure,
756 sources=sources,
757 skyInfo=skyInfo,
758 exposureId=idGenerator.catalog_id,
759 ccdInputs=ccdInputs,
760 sourceTableHandleDict=sourceTableHandleDict,
761 finalizedSourceTableHandleDict=finalizedSourceTableHandleDict,
762 finalVisitSummaryHandleDict=finalVisitSummaryHandleDict,
763 apCorrMap=apCorrMap,
764 )
765 # Strip HeavyFootprints to save space on disk
766 if self.config.doStripFootprints:
767 sources = outputs.outputSources
768 for source in sources[sources["parent"] != 0]:
769 source.setFootprint(None)
770 butlerQC.put(outputs, outputRefs)
772 def run(self, exposure, sources, skyInfo, exposureId, ccdInputs=None,
773 sourceTableHandleDict=None, finalizedSourceTableHandleDict=None, finalVisitSummaryHandleDict=None,
774 apCorrMap=None):
775 """Run measurement algorithms on the input exposure, and optionally populate the
776 resulting catalog with extra information.
778 Parameters
779 ----------
780 exposure : `lsst.afw.exposure.Exposure`
781 The input exposure on which measurements are to be performed.
782 sources : `lsst.afw.table.SourceCatalog`
783 A catalog built from the results of merged detections, or
784 deblender outputs.
785 parentCatalog : `lsst.afw.table.SourceCatalog`
786 Catalog of parent sources corresponding to sources.
787 skyInfo : `lsst.pipe.base.Struct`
788 A struct containing information about the position of the input exposure within
789 a `SkyMap`, the `SkyMap`, its `Wcs`, and its bounding box.
790 exposureId : `int` or `bytes`
791 Packed unique number or bytes unique to the input exposure.
792 ccdInputs : `lsst.afw.table.ExposureCatalog`, optional
793 Catalog containing information on the individual visits which went into making
794 the coadd.
795 sourceTableHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
796 Dict for sourceTable_visit handles (key is visit) for propagating flags.
797 These tables contain astrometry and photometry flags, and optionally PSF flags.
798 finalizedSourceTableHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
799 Dict for finalized_src_table handles (key is visit) for propagating flags.
800 These tables contain PSF flags from the finalized PSF estimation.
801 finalVisitSummaryHandleDict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`], optional
802 Dict for visit_summary handles (key is visit) for visit-level information.
803 These tables contain the WCS information of the single-visit input images.
804 apCorrMap : `lsst.afw.image.ApCorrMap`, optional
805 Aperture correction map attached to the ``exposure``. If None, it
806 will be read from the ``exposure``.
808 Returns
809 -------
810 results : `lsst.pipe.base.Struct`
811 Results of running measurement task. Will contain the catalog in the
812 sources attribute.
813 """
814 if self.config.doPropagateFlags:
815 # These mask planes may not be defined on the coadds always.
816 # We add the mask planes, which is a no-op if already defined.
817 for maskPlane in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
818 exposure.mask.addMaskPlane(maskPlane)
819 for maskPlane in self.config.measurement.plugins["base_PixelFlags"].masksFpCenter:
820 exposure.mask.addMaskPlane(maskPlane)
822 self.measurement.run(sources, exposure, exposureId=exposureId)
824 if self.config.doApCorr:
825 if apCorrMap is None:
826 apCorrMap = exposure.getInfo().getApCorrMap()
827 self.applyApCorr.run(
828 catalog=sources,
829 apCorrMap=apCorrMap,
830 )
832 # TODO DM-11568: this contiguous check-and-copy could go away if we
833 # reserve enough space during SourceDetection and/or SourceDeblend.
834 # NOTE: sourceSelectors require contiguous catalogs, so ensure
835 # contiguity now, so views are preserved from here on.
836 if not sources.isContiguous():
837 sources = sources.copy(deep=True)
839 if self.config.doRunCatalogCalculation:
840 self.catalogCalculation.run(sources)
842 self.setPrimaryFlags.run(sources, skyMap=skyInfo.skyMap, tractInfo=skyInfo.tractInfo,
843 patchInfo=skyInfo.patchInfo)
844 if self.config.doPropagateFlags:
845 self.propagateFlags.run(
846 sources,
847 ccdInputs,
848 sourceTableHandleDict,
849 finalizedSourceTableHandleDict,
850 finalVisitSummaryHandleDict,
851 )
853 results = Struct()
854 results.outputSources = sources
855 return results