Coverage for python/lsst/drp/tasks/assemble_coadd.py: 14%
670 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-11 02:25 -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/>.
22from __future__ import annotations
24__all__ = [
25 "AssembleCoaddTask",
26 "AssembleCoaddConnections",
27 "AssembleCoaddConfig",
28 "CompareWarpAssembleCoaddTask",
29 "CompareWarpAssembleCoaddConfig",
30]
32import copy
33import logging
34import warnings
36import lsstDebug
37import numpy
39import lsst.afw.geom as afwGeom
40import lsst.afw.image as afwImage
41import lsst.afw.math as afwMath
42import lsst.afw.table as afwTable
43import lsst.coadd.utils as coaddUtils
44import lsst.geom as geom
45import lsst.meas.algorithms as measAlg
46import lsst.pex.config as pexConfig
47import lsst.pex.exceptions as pexExceptions
48import lsst.pipe.base as pipeBase
49import lsst.utils as utils
50from lsst.meas.algorithms import AccumulatorMeanStack, MaskStreaksTask, ScaleVarianceTask, SourceDetectionTask
51from lsst.pipe.tasks.coaddBase import (
52 CoaddBaseTask,
53 makeSkyInfo,
54 removeMaskPlanes,
55 reorderAndPadList,
56 setRejectedMaskMapping,
57 subBBoxIter,
58)
59from lsst.pipe.tasks.healSparseMapping import HealSparseInputMapTask
60from lsst.pipe.tasks.interpImage import InterpImageTask
61from lsst.pipe.tasks.scaleZeroPoint import ScaleZeroPointTask
62from lsst.skymap import BaseSkyMap
63from lsst.utils.timer import timeMethod
65log = logging.getLogger(__name__)
68class AssembleCoaddConnections(
69 pipeBase.PipelineTaskConnections,
70 dimensions=("tract", "patch", "band", "skymap"),
71 defaultTemplates={
72 "inputCoaddName": "deep",
73 "outputCoaddName": "deep",
74 "warpType": "direct",
75 "warpTypeSuffix": "",
76 },
77):
78 inputWarps = pipeBase.connectionTypes.Input(
79 doc=(
80 "Input list of warps to be assemebled i.e. stacked."
81 "WarpType (e.g. direct, psfMatched) is controlled by the "
82 "warpType config parameter"
83 ),
84 name="{inputCoaddName}Coadd_{warpType}Warp",
85 storageClass="ExposureF",
86 dimensions=("tract", "patch", "skymap", "visit", "instrument"),
87 deferLoad=True,
88 multiple=True,
89 )
90 skyMap = pipeBase.connectionTypes.Input(
91 doc="Input definition of geometry/bbox and projection/wcs for coadded " "exposures",
92 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
93 storageClass="SkyMap",
94 dimensions=("skymap",),
95 )
96 selectedVisits = pipeBase.connectionTypes.Input(
97 doc="Selected visits to be coadded.",
98 name="{outputCoaddName}Visits",
99 storageClass="StructuredDataDict",
100 dimensions=("instrument", "tract", "patch", "skymap", "band"),
101 )
102 brightObjectMask = pipeBase.connectionTypes.PrerequisiteInput(
103 doc=(
104 "Input Bright Object Mask mask produced with external catalogs "
105 "to be applied to the mask plane BRIGHT_OBJECT."
106 ),
107 name="brightObjectMask",
108 storageClass="ObjectMaskCatalog",
109 dimensions=("tract", "patch", "skymap", "band"),
110 minimum=0,
111 )
112 coaddExposure = pipeBase.connectionTypes.Output(
113 doc="Output coadded exposure, produced by stacking input warps",
114 name="{outputCoaddName}Coadd{warpTypeSuffix}",
115 storageClass="ExposureF",
116 dimensions=("tract", "patch", "skymap", "band"),
117 )
118 nImage = pipeBase.connectionTypes.Output(
119 doc="Output image of number of input images per pixel",
120 name="{outputCoaddName}Coadd_nImage",
121 storageClass="ImageU",
122 dimensions=("tract", "patch", "skymap", "band"),
123 )
124 inputMap = pipeBase.connectionTypes.Output(
125 doc="Output healsparse map of input images",
126 name="{outputCoaddName}Coadd_inputMap",
127 storageClass="HealSparseMap",
128 dimensions=("tract", "patch", "skymap", "band"),
129 )
131 def __init__(self, *, config=None):
132 super().__init__(config=config)
134 if not config.doMaskBrightObjects:
135 self.prerequisiteInputs.remove("brightObjectMask")
137 if not config.doSelectVisits:
138 self.inputs.remove("selectedVisits")
140 if not config.doNImage:
141 self.outputs.remove("nImage")
143 if not self.config.doInputMap:
144 self.outputs.remove("inputMap")
146 if not self.config.doWriteArtifactMasks:
147 self.outputs.remove("artifactMasks")
150class AssembleCoaddConfig(
151 CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig, pipelineConnections=AssembleCoaddConnections
152):
153 warpType = pexConfig.ChoiceField(
154 doc="Warp name: one of 'direct' or 'psfMatched'",
155 dtype=str,
156 default="direct",
157 allowed={
158 "direct": "Weighted mean of directWarps, with outlier rejection",
159 "psfMatched": "Weighted mean of PSF-matched warps",
160 },
161 )
162 subregionSize = pexConfig.ListField(
163 dtype=int,
164 doc="Width, height of stack subregion size; "
165 "make small enough that a full stack of images will fit into memory "
166 " at once.",
167 length=2,
168 default=(2000, 2000),
169 )
170 statistic = pexConfig.Field(
171 dtype=str,
172 doc="Main stacking statistic for aggregating over the epochs.",
173 default="MEANCLIP",
174 )
175 doOnlineForMean = pexConfig.Field(
176 dtype=bool,
177 doc='Perform online coaddition when statistic="MEAN" to save memory?',
178 default=False,
179 )
180 sigmaClip = pexConfig.Field(
181 dtype=float,
182 doc="Sigma for outlier rejection; ignored if non-clipping statistic " "selected.",
183 default=3.0,
184 )
185 clipIter = pexConfig.Field(
186 dtype=int,
187 doc="Number of iterations of outlier rejection; ignored if " "non-clipping statistic selected.",
188 default=2,
189 )
190 calcErrorFromInputVariance = pexConfig.Field(
191 dtype=bool,
192 doc="Calculate coadd variance from input variance by stacking "
193 "statistic. Passed to "
194 "StatisticsControl.setCalcErrorFromInputVariance()",
195 default=True,
196 )
197 doScaleZeroPoint = pexConfig.Field(
198 dtype=bool,
199 doc="Scale the photometric zero point of the coadd temp exposures "
200 "such that the magnitude zero point results in a flux in nJy.",
201 default=False,
202 deprecated="Now that visits are scaled to nJy it is no longer necessary or "
203 "recommended to scale the zero point, so this will be removed "
204 "after v29.",
205 )
206 scaleZeroPoint = pexConfig.ConfigurableField(
207 target=ScaleZeroPointTask,
208 doc="Task to adjust the photometric zero point of the coadd temp " "exposures",
209 deprecated="Now that visits are scaled to nJy it is no longer necessary or "
210 "recommended to scale the zero point, so this will be removed "
211 "after v29.",
212 )
213 doInterp = pexConfig.Field(
214 doc="Interpolate over NaN pixels? Also extrapolate, if necessary, but " "the results are ugly.",
215 dtype=bool,
216 default=True,
217 )
218 interpImage = pexConfig.ConfigurableField(
219 target=InterpImageTask,
220 doc="Task to interpolate (and extrapolate) over NaN pixels",
221 )
222 doWrite = pexConfig.Field(
223 doc="Persist coadd?",
224 dtype=bool,
225 default=True,
226 )
227 doWriteArtifactMasks = pexConfig.Field(
228 doc="Persist artifact masks? Should be True for CompareWarp only.",
229 dtype=bool,
230 default=False,
231 )
232 doNImage = pexConfig.Field(
233 doc="Create image of number of contributing exposures for each pixel",
234 dtype=bool,
235 default=False,
236 )
237 maskPropagationThresholds = pexConfig.DictField(
238 keytype=str,
239 itemtype=float,
240 doc=(
241 "Threshold (in fractional weight) of rejection at which we "
242 "propagate a mask plane to the coadd; that is, we set the mask "
243 "bit on the coadd if the fraction the rejected frames "
244 "would have contributed exceeds this value."
245 ),
246 default={"SAT": 0.1},
247 )
248 removeMaskPlanes = pexConfig.ListField(
249 dtype=str,
250 default=["NOT_DEBLENDED"],
251 doc="Mask planes to remove before coadding",
252 )
253 doMaskBrightObjects = pexConfig.Field(
254 dtype=bool,
255 default=False,
256 doc="Set mask and flag bits for bright objects?",
257 )
258 brightObjectMaskName = pexConfig.Field(
259 dtype=str,
260 default="BRIGHT_OBJECT",
261 doc="Name of mask bit used for bright objects",
262 )
263 coaddPsf = pexConfig.ConfigField(
264 doc="Configuration for CoaddPsf",
265 dtype=measAlg.CoaddPsfConfig,
266 )
267 doAttachTransmissionCurve = pexConfig.Field(
268 dtype=bool,
269 default=False,
270 optional=False,
271 doc=(
272 "Attach a piecewise TransmissionCurve for the coadd? "
273 "(requires all input Exposures to have TransmissionCurves)."
274 ),
275 )
276 hasFakes = pexConfig.Field(
277 dtype=bool,
278 default=False,
279 doc="Should be set to True if fake sources have been inserted into the input data.",
280 )
281 doSelectVisits = pexConfig.Field(
282 doc="Coadd only visits selected by a SelectVisitsTask",
283 dtype=bool,
284 default=False,
285 )
286 doInputMap = pexConfig.Field(
287 doc="Create a bitwise map of coadd inputs",
288 dtype=bool,
289 default=False,
290 )
291 inputMapper = pexConfig.ConfigurableField(
292 doc="Input map creation subtask.",
293 target=HealSparseInputMapTask,
294 )
296 def setDefaults(self):
297 super().setDefaults()
298 self.badMaskPlanes = ["NO_DATA", "BAD", "SAT", "EDGE"]
299 self.coaddPsf.cacheSize = 0
301 def validate(self):
302 super().validate()
303 if self.doInterp and self.statistic not in ["MEAN", "MEDIAN", "MEANCLIP", "VARIANCE", "VARIANCECLIP"]:
304 raise pexConfig.FieldValidationError(
305 self.__class__.doInterp,
306 self,
307 f"Must set doInterp=False for statistic={self.statistic}, which does not "
308 "compute and set a non-zero coadd variance estimate.",
309 )
311 unstackableStats = ["NOTHING", "ERROR", "ORMASK"]
312 if not hasattr(afwMath.Property, self.statistic) or self.statistic in unstackableStats:
313 stackableStats = [
314 str(k) for k in afwMath.Property.__members__.keys() if str(k) not in unstackableStats
315 ]
316 raise pexConfig.FieldValidationError(
317 self.__class__.statistic,
318 self,
319 f"statistic {self.statistic} is not allowed. Please choose one of {stackableStats}.",
320 )
322 # Admittedly, it's odd for a parent class to condition on a child class
323 # but such is the case until the CompareWarp refactor in DM-38630.
324 if self.doWriteArtifactMasks and not isinstance(self, CompareWarpAssembleCoaddConfig):
325 raise pexConfig.FieldValidationError(
326 self.__class__.doWriteArtifactMasks,
327 self,
328 "doWriteArtifactMasks is only valid for CompareWarpAssembleCoaddConfig.",
329 )
332class AssembleCoaddTask(CoaddBaseTask, pipeBase.PipelineTask):
333 """Assemble a coadded image from a set of warps.
335 Each Warp that goes into a coadd will have its flux calibrated to
336 nJy. WarpType may be one of 'direct' or
337 'psfMatched', and the boolean configs `config.makeDirect` and
338 `config.makePsfMatched` set which of the warp types will be coadded.
339 The coadd is computed as a mean with optional outlier rejection.
340 Criteria for outlier rejection are set in `AssembleCoaddConfig`.
341 Finally, Warps can have bad 'NaN' pixels which received no input from the
342 source calExps. We interpolate over these bad (NaN) pixels.
344 `AssembleCoaddTask` uses several sub-tasks. These are
346 - `~lsst.pipe.tasks.ScaleZeroPointTask`
347 - create and use an ``imageScaler`` object to scale the photometric
348 zeropoint for each Warp (deprecated and will be removed in DM-49083).
349 - `~lsst.pipe.tasks.InterpImageTask`
350 - interpolate across bad pixels (NaN) in the final coadd
352 You can retarget these subtasks if you wish.
354 Raises
355 ------
356 RuntimeError
357 Raised if unable to define mask plane for bright objects.
359 Notes
360 -----
361 Debugging:
362 `AssembleCoaddTask` has no debug variables of its own. Some of the
363 subtasks may support `~lsst.base.lsstDebug` variables. See the
364 documentation for the subtasks for further information.
365 """
367 ConfigClass = AssembleCoaddConfig
368 _DefaultName = "assembleCoadd"
370 _doUsePsfMatchedPolygons: bool = False
371 """Use ValidPolygons from shrunk Psf-Matched Calexps?
373 This needs to be set to True by child classes that use compare Psf-Matched
374 warps to non Psf-Matched warps.
375 """
377 def __init__(self, *args, **kwargs):
378 # TODO: DM-17415 better way to handle previously allowed passed args
379 # e.g.`AssembleCoaddTask(config)`
380 if args:
381 argNames = ["config", "name", "parentTask", "log"]
382 kwargs.update({k: v for k, v in zip(argNames, args)})
383 warnings.warn(
384 "AssembleCoadd received positional args, and casting them as kwargs: %s. "
385 "PipelineTask will not take positional args" % argNames,
386 FutureWarning,
387 stacklevel=2,
388 )
390 super().__init__(**kwargs)
391 self.makeSubtask("interpImage")
392 if self.config.doScaleZeroPoint:
393 # Remove completely in DM-49083
394 self.makeSubtask("scaleZeroPoint")
396 if self.config.doMaskBrightObjects:
397 mask = afwImage.Mask()
398 try:
399 self.brightObjectBitmask = 1 << mask.addMaskPlane(self.config.brightObjectMaskName)
400 except pexExceptions.LsstCppException:
401 raise RuntimeError(
402 "Unable to define mask plane for bright objects; planes used are %s"
403 % mask.getMaskPlaneDict().keys()
404 )
405 del mask
407 if self.config.doInputMap:
408 self.makeSubtask("inputMapper")
410 self.warpType = self.config.warpType
412 def runQuantum(self, butlerQC, inputRefs, outputRefs):
413 inputData = butlerQC.get(inputRefs)
415 # Construct skyInfo expected by run
416 # Do not remove skyMap from inputData in case _makeSupplementaryData
417 # needs it
418 skyMap = inputData["skyMap"]
419 outputDataId = butlerQC.quantum.dataId
421 inputData["skyInfo"] = makeSkyInfo(
422 skyMap, tractId=outputDataId["tract"], patchId=outputDataId["patch"]
423 )
425 if self.config.doSelectVisits:
426 warpRefList = self.filterWarps(inputData["inputWarps"], inputData["selectedVisits"])
427 else:
428 warpRefList = inputData["inputWarps"]
430 # Although psfMatchedWarps are specifically required by only by a
431 # specific set of a subclass, and since runQuantum is not overridden,
432 # the selection and filtering must happen here on a conditional basis.
433 # Otherwise, the elements in the various lists will not line up.
434 # This design cries out for a refactor, which is planned in DM-38630.
435 if self._doUsePsfMatchedPolygons:
436 if self.config.doSelectVisits:
437 psfMatchedWarpRefList = self.filterWarps(
438 inputData["psfMatchedWarps"],
439 inputData["selectedVisits"],
440 )
441 else:
442 psfMatchedWarpRefList = inputData["psfMatchedWarps"]
443 else:
444 psfMatchedWarpRefList = []
446 inputs = self.prepareInputs(warpRefList, inputData["skyInfo"].bbox, psfMatchedWarpRefList)
447 self.log.info("Found %d %s", len(inputs.warpRefList), self.getTempExpDatasetName(self.warpType))
448 if len(inputs.warpRefList) == 0:
449 raise pipeBase.NoWorkFound("No coadd temporary exposures found")
451 supplementaryData = self._makeSupplementaryData(butlerQC, inputRefs, outputRefs)
452 retStruct = self.run(
453 inputData["skyInfo"],
454 warpRefList=inputs.warpRefList,
455 imageScalerList=inputs.imageScalerList,
456 weightList=inputs.weightList,
457 psfMatchedWarpRefList=inputs.psfMatchedWarpRefList,
458 supplementaryData=supplementaryData,
459 )
461 inputData.setdefault("brightObjectMask", None)
462 if self.config.doMaskBrightObjects and inputData["brightObjectMask"] is None:
463 log.warning("doMaskBrightObjects is set to True, but brightObjectMask not loaded")
464 self.processResults(retStruct.coaddExposure, inputData["brightObjectMask"], outputDataId)
466 if self.config.doWriteArtifactMasks:
467 artifactMasksRefList = reorderAndPadList(
468 outputRefs.artifactMasks,
469 [ref.dataId for ref in outputRefs.artifactMasks],
470 [ref.dataId for ref in inputs.warpRefList],
471 )
472 for altMask, warpRef, outputRef in zip(
473 retStruct.altMaskList, inputs.warpRefList, artifactMasksRefList, strict=True
474 ):
475 mask = warpRef.get(component="mask", parameters={"bbox": retStruct.coaddExposure.getBBox()})
476 self.applyAltMaskPlanes(mask, altMask)
477 butlerQC.put(mask, outputRef)
479 if self.config.doWrite:
480 butlerQC.put(retStruct, outputRefs)
482 return retStruct
484 def processResults(self, coaddExposure, brightObjectMasks=None, dataId=None):
485 """Interpolate over missing data and mask bright stars.
487 Parameters
488 ----------
489 coaddExposure : `lsst.afw.image.Exposure`
490 The coadded exposure to process.
491 brightObjectMasks : `lsst.afw.table` or `None`, optional
492 Table of bright objects to mask.
493 dataId : `lsst.daf.butler.DataId` or `None`, optional
494 Data identification.
495 """
496 if self.config.doInterp:
497 self.interpImage.run(coaddExposure.getMaskedImage(), planeName="NO_DATA")
498 # The variance must be positive; work around for DM-3201.
499 varArray = coaddExposure.variance.array
500 with numpy.errstate(invalid="ignore"):
501 varArray[:] = numpy.where(varArray > 0, varArray, numpy.inf)
503 if self.config.doMaskBrightObjects:
504 self.setBrightObjectMasks(coaddExposure, brightObjectMasks, dataId)
506 def _makeSupplementaryData(self, butlerQC, inputRefs, outputRefs):
507 """Make additional inputs to run() specific to subclasses (Gen3).
509 Duplicates interface of `runQuantum` method.
510 Available to be implemented by subclasses only if they need the
511 coadd dataRef/handle for performing preliminary processing before
512 assembling the coadd.
514 Parameters
515 ----------
516 butlerQC : `~lsst.pipe.base.QuantumContext`
517 Gen3 Butler object for fetching additional data products before
518 running the Task specialized for quantum being processed.
519 inputRefs : `~lsst.pipe.base.InputQuantizedConnection`
520 Attributes are the names of the connections describing input
521 dataset types. Values are DatasetRefs that task consumes for the
522 corresponding dataset type. DataIds are guaranteed to match data
523 objects in ``inputData``.
524 outputRefs : `~lsst.pipe.base.OutputQuantizedConnection`
525 Attributes are the names of the connections describing output
526 dataset types. Values are DatasetRefs that task is to produce
527 for the corresponding dataset type.
528 """
529 return pipeBase.Struct()
531 def prepareInputs(self, refList, coadd_bbox, psfMatchedWarpRefList=None):
532 """Prepare the input warps for coaddition by measuring the weight for
533 each warp.
535 Before coadding these Warps together compute the weight for each
536 Warp.
538 Parameters
539 ----------
540 refList : `list`
541 List of dataset handles (data references) to warp.
542 psfMatchedWarpRefList : `list` | None, optional
543 List of dataset handles (data references) to psfMatchedWarp.
545 Returns
546 -------
547 result : `~lsst.pipe.base.Struct`
548 Results as a struct with attributes:
550 ``warpRefList``
551 `list` of dataset handles (data references) to warp.
552 ``weightList``
553 `list` of weightings.
554 ``imageScalerList``
555 `list` of image scalers.
556 Deprecated and will be removed in DM-49083.
557 """
558 statsCtrl = afwMath.StatisticsControl()
559 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
560 statsCtrl.setNumIter(self.config.clipIter)
561 statsCtrl.setAndMask(self.getBadPixelMask())
562 statsCtrl.setNanSafe(True)
563 # compute warpRefList: a list of warpRef that actually exist
564 # and weightList: a list of the weight of the associated coadd warp
565 # and imageScalerList: a list of scale factors for the associated coadd
566 # warp. (deprecated and will be removed in DM-49083)
567 warpRefList = []
568 weightList = []
569 # Remove in DM-49083
570 imageScalerList = []
571 outputPsfMatchedWarpRefList = []
573 # Convert psfMatchedWarpRefList to a dict, so we can look them up by
574 # their dataId.
575 # Note: Some warps may not have a corresponding psfMatchedWarp, which
576 # could have happened due to a failure in the PSF matching, or not
577 # having enough good pixels to support the PSF-matching kernel.
579 if psfMatchedWarpRefList is None:
580 psfMatchedWarpRefDict = {ref.dataId: None for ref in refList}
581 else:
582 psfMatchedWarpRefDict = {ref.dataId: ref for ref in psfMatchedWarpRefList}
584 warpName = self.getTempExpDatasetName(self.warpType)
585 for warpRef in refList:
586 warp = warpRef.get(parameters={"bbox": coadd_bbox})
587 # Ignore any input warp that is empty of data
588 if numpy.isnan(warp.image.array).all():
589 continue
590 maskedImage = warp.getMaskedImage()
592 # Allow users to scale the zero point for backward
593 # compatibility. Remove imageScalarList in DM-49083.
594 if self.config.doScaleZeroPoint:
595 imageScaler = self.scaleZeroPoint.computeImageScaler(
596 exposure=warp,
597 dataRef=warpRef, # FIXME
598 )
599 try:
600 imageScaler.scaleMaskedImage(maskedImage)
601 except Exception as e:
602 self.log.warning("Scaling failed for %s (skipping it): %s", warpRef.dataId, e)
603 continue
604 imageScalerList.append(imageScaler)
605 else:
606 imageScalerList.append(None)
607 if "BUNIT" not in warp.metadata:
608 raise ValueError(f"Warp {warpRef.dataId} has no BUNIT metadata")
609 if warp.metadata["BUNIT"] != "nJy":
610 raise ValueError(
611 f"Warp {warpRef.dataId} has BUNIT {warp.metadata['BUNIT']}, expected nJy"
612 )
613 statObj = afwMath.makeStatistics(
614 maskedImage.getVariance(), maskedImage.getMask(), afwMath.MEANCLIP, statsCtrl
615 )
616 meanVar, meanVarErr = statObj.getResult(afwMath.MEANCLIP)
617 if meanVar > 0.0 and numpy.isfinite(meanVar):
618 weight = 1.0 / float(meanVar)
619 self.log.info("Weight of %s %s = %0.3f", warpName, warpRef.dataId, weight)
620 else:
621 self.log.warning("Non-finite weight for %s: skipping", warpRef.dataId)
622 continue
624 del maskedImage
625 del warp
627 warpRefList.append(warpRef)
628 weightList.append(weight)
630 dataId = warpRef.dataId
631 psfMatchedWarpRef = psfMatchedWarpRefDict.get(dataId, None)
632 outputPsfMatchedWarpRefList.append(psfMatchedWarpRef)
634 return pipeBase.Struct(
635 warpRefList=warpRefList,
636 weightList=weightList,
637 imageScalerList=imageScalerList,
638 psfMatchedWarpRefList=outputPsfMatchedWarpRefList,
639 )
641 def prepareStats(self, mask=None):
642 """Prepare the statistics for coadding images.
644 Parameters
645 ----------
646 mask : `int`, optional
647 Bit mask value to exclude from coaddition.
649 Returns
650 -------
651 stats : `~lsst.pipe.base.Struct`
652 Statistics as a struct with attributes:
654 ``statsCtrl``
655 Statistics control object for coadd
656 (`~lsst.afw.math.StatisticsControl`).
657 ``statsFlags``
658 Statistic for coadd (`~lsst.afw.math.Property`).
659 """
660 if mask is None:
661 mask = self.getBadPixelMask()
662 statsCtrl = afwMath.StatisticsControl()
663 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
664 statsCtrl.setNumIter(self.config.clipIter)
665 statsCtrl.setAndMask(mask)
666 statsCtrl.setNanSafe(True)
667 statsCtrl.setWeighted(True)
668 statsCtrl.setCalcErrorFromInputVariance(self.config.calcErrorFromInputVariance)
669 for plane, threshold in self.config.maskPropagationThresholds.items():
670 bit = afwImage.Mask.getMaskPlane(plane)
671 statsCtrl.setMaskPropagationThreshold(bit, threshold)
672 statsFlags = afwMath.stringToStatisticsProperty(self.config.statistic)
673 return pipeBase.Struct(ctrl=statsCtrl, flags=statsFlags)
675 @timeMethod
676 def run(
677 self,
678 skyInfo,
679 *,
680 warpRefList,
681 imageScalerList,
682 weightList,
683 psfMatchedWarpRefList=None,
684 altMaskList=None,
685 mask=None,
686 supplementaryData=None,
687 ):
688 """Assemble a coadd from input warps.
690 Assemble the coadd using the provided list of coaddTempExps. Since
691 the full coadd covers a patch (a large area), the assembly is
692 performed over small areas on the image at a time in order to
693 conserve memory usage. Iterate over subregions within the outer
694 bbox of the patch using `assembleSubregion` to stack the corresponding
695 subregions from the coaddTempExps with the statistic specified.
696 Set the edge bits the coadd mask based on the weight map.
698 Parameters
699 ----------
700 skyInfo : `~lsst.pipe.base.Struct`
701 Struct with geometric information about the patch.
702 warpRefList : `list`
703 List of dataset handles (data references) to Warps
704 (previously called CoaddTempExps).
705 imageScalerList : `list`
706 List of image scalers. Deprecated and will be removed after v29
707 in DM-49083.
708 weightList : `list`
709 List of weights.
710 psfMatchedWarpRefList : `list`, optional
711 List of dataset handles (data references) to psfMatchedWarps.
712 altMaskList : `list`, optional
713 List of alternate masks to use rather than those stored with
714 warp.
715 mask : `int`, optional
716 Bit mask value to exclude from coaddition.
717 supplementaryData : `~lsst.pipe.base.Struct`, optional
718 Struct with additional data products needed to assemble coadd.
719 Only used by subclasses that implement ``_makeSupplementaryData``
720 and override `run`.
722 Returns
723 -------
724 result : `~lsst.pipe.base.Struct`
725 Results as a struct with attributes:
727 ``coaddExposure``
728 Coadded exposure (`~lsst.afw.image.Exposure`).
729 ``nImage``
730 Exposure count image (`~lsst.afw.image.Image`), if requested.
731 ``inputMap``
732 Bit-wise map of inputs, if requested.
733 ``warpRefList``
734 Input list of dataset handles (data refs) to the warps
735 (`~lsst.daf.butler.DeferredDatasetHandle`) (unmodified).
736 ``imageScalerList``
737 Input list of image scalers (`list`) (unmodified).
738 Deprecated and will be removed after v29 in DM-49083.
739 ``weightList``
740 Input list of weights (`list`) (unmodified).
742 Raises
743 ------
744 lsst.pipe.base.NoWorkFound
745 Raised if no dataset handles (data references) are provided.
746 """
747 warpName = self.getTempExpDatasetName(self.warpType)
748 self.log.info("Assembling %s %s", len(warpRefList), warpName)
749 if not warpRefList:
750 raise pipeBase.NoWorkFound("No exposures provided for co-addition.")
752 stats = self.prepareStats(mask=mask)
754 if altMaskList is None:
755 altMaskList = [None] * len(warpRefList)
757 coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs)
758 # Deprecated, keep only the `else` branch in DM-49083.
759 if self.config.doScaleZeroPoint:
760 coaddExposure.setPhotoCalib(self.scaleZeroPoint.getPhotoCalib())
761 else:
762 coaddExposure.setPhotoCalib(afwImage.PhotoCalib(1.0))
763 coaddExposure.getInfo().setCoaddInputs(self.inputRecorder.makeCoaddInputs())
764 self.assembleMetadata(coaddExposure, warpRefList, weightList, psfMatchedWarpRefList)
765 coaddMaskedImage = coaddExposure.getMaskedImage()
766 subregionSizeArr = self.config.subregionSize
767 subregionSize = geom.Extent2I(subregionSizeArr[0], subregionSizeArr[1])
768 # if nImage is requested, create a zero one which can be passed to
769 # assembleSubregion.
770 if self.config.doNImage:
771 nImage = afwImage.ImageU(skyInfo.bbox)
772 else:
773 nImage = None
774 # If inputMap is requested, create the initial version that can be
775 # masked in assembleSubregion.
776 if self.config.doInputMap:
777 self.inputMapper.build_ccd_input_map(
778 skyInfo.bbox, skyInfo.wcs, coaddExposure.getInfo().getCoaddInputs().ccds
779 )
781 if self.config.doOnlineForMean and self.config.statistic == "MEAN":
782 try:
783 self.assembleOnlineMeanCoadd(
784 coaddExposure,
785 warpRefList,
786 imageScalerList,
787 weightList,
788 altMaskList,
789 stats.ctrl,
790 nImage=nImage,
791 )
792 except Exception as e:
793 self.log.exception("Cannot compute online coadd %s", e)
794 raise
795 else:
796 for subBBox in subBBoxIter(skyInfo.bbox, subregionSize):
797 try:
798 self.assembleSubregion(
799 coaddExposure,
800 subBBox,
801 warpRefList,
802 imageScalerList,
803 weightList,
804 altMaskList,
805 stats.flags,
806 stats.ctrl,
807 nImage=nImage,
808 )
809 except Exception as e:
810 self.log.exception("Cannot compute coadd %s: %s", subBBox, e)
811 raise
813 # If inputMap is requested, we must finalize the map after the
814 # accumulation.
815 if self.config.doInputMap:
816 self.inputMapper.finalize_ccd_input_map_mask()
817 inputMap = self.inputMapper.ccd_input_map
818 else:
819 inputMap = None
821 self.setInexactPsf(coaddMaskedImage.getMask())
822 # Despite the name, the following doesn't really deal with "EDGE"
823 # pixels: it identifies pixels that didn't receive any unmasked inputs
824 # (as occurs around the edge of the field).
825 coaddUtils.setCoaddEdgeBits(coaddMaskedImage.getMask(), coaddMaskedImage.getVariance())
826 return pipeBase.Struct(
827 coaddExposure=coaddExposure,
828 nImage=nImage,
829 warpRefList=warpRefList,
830 imageScalerList=imageScalerList,
831 weightList=weightList,
832 inputMap=inputMap,
833 altMaskList=altMaskList,
834 )
836 def assembleMetadata(self, coaddExposure, warpRefList, weightList, psfMatchedWarpRefList=None):
837 """Set the metadata for the coadd.
839 This basic implementation sets the filter from the first input.
841 Parameters
842 ----------
843 coaddExposure : `lsst.afw.image.Exposure`
844 The target exposure for the coadd.
845 warpRefList : `list`
846 List of dataset handles (data references) to warp.
847 weightList : `list`
848 List of weights.
849 psfMatchedWarpRefList : `list` | None, optional
850 List of dataset handles (data references) to psfMatchedWarps.
852 Raises
853 ------
854 AssertionError
855 Raised if there is a length mismatch.
856 """
857 assert len(warpRefList) == len(weightList), "Length mismatch"
859 if psfMatchedWarpRefList:
860 assert len(warpRefList) == len(psfMatchedWarpRefList), "Length mismatch"
862 # We load a single pixel of each coaddTempExp, because we just want to
863 # get at the metadata (and we need more than just the PropertySet that
864 # contains the header), which is not possible with the current butler
865 # (see #2777).
866 bbox = geom.Box2I(coaddExposure.getBBox().getMin(), geom.Extent2I(1, 1))
868 warpList = [warpRef.get(parameters={"bbox": bbox}) for warpRef in warpRefList]
870 numCcds = sum(len(warp.getInfo().getCoaddInputs().ccds) for warp in warpList)
872 # Set the coadd FilterLabel to the band of the first input exposure:
873 # Coadds are calibrated, so the physical label is now meaningless.
874 coaddExposure.setFilter(afwImage.FilterLabel(warpList[0].getFilter().bandLabel))
875 coaddInputs = coaddExposure.getInfo().getCoaddInputs()
876 coaddInputs.ccds.reserve(numCcds)
877 coaddInputs.visits.reserve(len(warpList))
879 # Set the exposure units to nJy
880 # No need to check after DM-49083.
881 if not self.config.doScaleZeroPoint:
882 coaddExposure.metadata["BUNIT"] = "nJy"
884 # psfMatchedWarpRefList should be empty except in CompareWarpCoadd.
885 if self._doUsePsfMatchedPolygons:
886 # Set validPolygons for warp before addVisitToCoadd
887 self._setValidPolygons(warpList, psfMatchedWarpRefList)
889 for warp, weight in zip(warpList, weightList):
890 self.inputRecorder.addVisitToCoadd(coaddInputs, warp, weight)
892 coaddInputs.visits.sort()
893 coaddInputs.ccds.sort()
894 if self.warpType == "psfMatched":
895 # The modelPsf BBox for a psfMatchedWarp/coaddTempExp was
896 # dynamically defined by ModelPsfMatchTask as the square box
897 # bounding its spatially-variable, pre-matched WarpedPsf.
898 # Likewise, set the PSF of a PSF-Matched Coadd to the modelPsf
899 # having the maximum width (sufficient because square)
900 modelPsfList = [warp.getPsf() for warp in warpList]
901 modelPsfWidthList = [
902 modelPsf.computeBBox(modelPsf.getAveragePosition()).getWidth() for modelPsf in modelPsfList
903 ]
904 psf = modelPsfList[modelPsfWidthList.index(max(modelPsfWidthList))]
905 else:
906 psf = measAlg.CoaddPsf(
907 coaddInputs.ccds, coaddExposure.getWcs(), self.config.coaddPsf.makeControl()
908 )
909 coaddExposure.setPsf(psf)
910 apCorrMap = measAlg.makeCoaddApCorrMap(
911 coaddInputs.ccds, coaddExposure.getBBox(afwImage.PARENT), coaddExposure.getWcs()
912 )
913 coaddExposure.getInfo().setApCorrMap(apCorrMap)
914 if self.config.doAttachTransmissionCurve:
915 transmissionCurve = measAlg.makeCoaddTransmissionCurve(coaddExposure.getWcs(), coaddInputs.ccds)
916 coaddExposure.getInfo().setTransmissionCurve(transmissionCurve)
918 def assembleSubregion(
919 self,
920 coaddExposure,
921 bbox,
922 warpRefList,
923 imageScalerList,
924 weightList,
925 altMaskList,
926 statsFlags,
927 statsCtrl,
928 nImage=None,
929 ):
930 """Assemble the coadd for a sub-region.
932 For each coaddTempExp, check for (and swap in) an alternative mask
933 if one is passed. Remove mask planes listed in
934 `config.removeMaskPlanes`. Finally, stack the actual exposures using
935 `lsst.afw.math.statisticsStack` with the statistic specified by
936 statsFlags. Typically, the statsFlag will be one of lsst.afw.math.MEAN
937 for a mean-stack or `lsst.afw.math.MEANCLIP` for outlier rejection
938 using an N-sigma clipped mean where N and iterations are specified by
939 statsCtrl. Assign the stacked subregion back to the coadd.
941 Parameters
942 ----------
943 coaddExposure : `lsst.afw.image.Exposure`
944 The target exposure for the coadd.
945 bbox : `lsst.geom.Box`
946 Sub-region to coadd.
947 warpRefList : `list`
948 List of dataset handles (data references) to warp.
949 imageScalerList : `list`
950 List of image scalers.
951 Deprecated and will be removed after v29 in DM-49083.
952 weightList : `list`
953 List of weights.
954 altMaskList : `list`
955 List of alternate masks to use rather than those stored with
956 warp, or None. Each element is dict with keys = mask plane
957 name to which to add the spans.
958 statsFlags : `lsst.afw.math.Property`
959 Property object for statistic for coadd.
960 statsCtrl : `lsst.afw.math.StatisticsControl`
961 Statistics control object for coadd.
962 nImage : `lsst.afw.image.ImageU`, optional
963 Keeps track of exposure count for each pixel.
964 """
965 self.log.info("Stacking subregion %s", bbox)
967 coaddExposure.mask.addMaskPlane("REJECTED")
968 coaddExposure.mask.addMaskPlane("CLIPPED")
969 coaddExposure.mask.addMaskPlane("SENSOR_EDGE")
970 maskMap = setRejectedMaskMapping(statsCtrl)
971 clipped = afwImage.Mask.getPlaneBitMask("CLIPPED")
972 maskedImageList = []
973 if nImage is not None:
974 subNImage = afwImage.ImageU(bbox.getWidth(), bbox.getHeight())
975 for warpRef, imageScaler, altMask in zip(warpRefList, imageScalerList, altMaskList):
976 exposure = warpRef.get(parameters={"bbox": bbox})
978 maskedImage = exposure.getMaskedImage()
979 mask = maskedImage.getMask()
980 if altMask is not None:
981 self.applyAltMaskPlanes(mask, altMask)
982 # Remove in DM-49083
983 if imageScaler is not None:
984 imageScaler.scaleMaskedImage(maskedImage)
986 # Add 1 for each pixel which is not excluded by the exclude mask.
987 # In legacyCoadd, pixels may also be excluded by
988 # afwMath.statisticsStack.
989 if nImage is not None:
990 subNImage.getArray()[maskedImage.getMask().getArray() & statsCtrl.getAndMask() == 0] += 1
991 if self.config.removeMaskPlanes:
992 removeMaskPlanes(maskedImage.mask, self.config.removeMaskPlanes, logger=self.log)
993 maskedImageList.append(maskedImage)
995 if self.config.doInputMap:
996 visit = exposure.getInfo().getCoaddInputs().visits[0].getId()
997 self.inputMapper.mask_warp_bbox(bbox, visit, mask, statsCtrl.getAndMask())
999 with self.timer("stack"):
1000 coaddSubregion = afwMath.statisticsStack(
1001 maskedImageList,
1002 statsFlags,
1003 statsCtrl,
1004 weightList,
1005 clipped, # also set output to CLIPPED if sigma-clipped
1006 maskMap,
1007 )
1008 coaddExposure.maskedImage.assign(coaddSubregion, bbox)
1009 if nImage is not None:
1010 nImage.assign(subNImage, bbox)
1012 def assembleOnlineMeanCoadd(
1013 self, coaddExposure, warpRefList, imageScalerList, weightList, altMaskList, statsCtrl, nImage=None
1014 ):
1015 """Assemble the coadd using the "online" method.
1017 This method takes a running sum of images and weights to save memory.
1018 It only works for MEAN statistics.
1020 Parameters
1021 ----------
1022 coaddExposure : `lsst.afw.image.Exposure`
1023 The target exposure for the coadd.
1024 warpRefList : `list`
1025 List of dataset handles (data references) to warp.
1026 imageScalerList : `list`
1027 List of image scalers.
1028 Deprecated and will be removed after v29 in DM-49083.
1029 weightList : `list`
1030 List of weights.
1031 altMaskList : `list`
1032 List of alternate masks to use rather than those stored with
1033 warp, or None. Each element is dict with keys = mask plane
1034 name to which to add the spans.
1035 statsCtrl : `lsst.afw.math.StatisticsControl`
1036 Statistics control object for coadd.
1037 nImage : `lsst.afw.image.ImageU`, optional
1038 Keeps track of exposure count for each pixel.
1039 """
1040 self.log.debug("Computing online coadd.")
1042 coaddExposure.mask.addMaskPlane("REJECTED")
1043 coaddExposure.mask.addMaskPlane("CLIPPED")
1044 coaddExposure.mask.addMaskPlane("SENSOR_EDGE")
1045 maskMap = setRejectedMaskMapping(statsCtrl)
1046 thresholdDict = AccumulatorMeanStack.stats_ctrl_to_threshold_dict(statsCtrl)
1048 bbox = coaddExposure.maskedImage.getBBox()
1050 stacker = AccumulatorMeanStack(
1051 coaddExposure.image.array.shape,
1052 statsCtrl.getAndMask(),
1053 mask_threshold_dict=thresholdDict,
1054 mask_map=maskMap,
1055 no_good_pixels_mask=statsCtrl.getNoGoodPixelsMask(),
1056 calc_error_from_input_variance=self.config.calcErrorFromInputVariance,
1057 compute_n_image=(nImage is not None),
1058 )
1060 for warpRef, imageScaler, altMask, weight in zip(
1061 warpRefList, imageScalerList, altMaskList, weightList
1062 ):
1063 exposure = warpRef.get(parameters={"bbox": bbox})
1064 maskedImage = exposure.getMaskedImage()
1065 mask = maskedImage.getMask()
1066 if altMask is not None:
1067 self.applyAltMaskPlanes(mask, altMask)
1068 # Remove in DM-49083.
1069 if imageScaler is not None:
1070 imageScaler.scaleMaskedImage(maskedImage)
1071 if self.config.removeMaskPlanes:
1072 removeMaskPlanes(maskedImage.mask, self.config.removeMaskPlanes, logger=self.log)
1074 stacker.add_masked_image(maskedImage, weight=weight)
1076 if self.config.doInputMap:
1077 visit = exposure.getInfo().getCoaddInputs().visits[0].getId()
1078 self.inputMapper.mask_warp_bbox(bbox, visit, mask, statsCtrl.getAndMask())
1080 stacker.fill_stacked_masked_image(coaddExposure.maskedImage)
1082 if nImage is not None:
1083 nImage.array[:, :] = stacker.n_image
1085 def applyAltMaskPlanes(self, mask, altMaskSpans):
1086 """Apply in place alt mask formatted as SpanSets to a mask.
1088 Parameters
1089 ----------
1090 mask : `lsst.afw.image.Mask`
1091 Original mask.
1092 altMaskSpans : `dict`
1093 SpanSet lists to apply. Each element contains the new mask
1094 plane name (e.g. "CLIPPED and/or "NO_DATA") as the key,
1095 and list of SpanSets to apply to the mask.
1097 Returns
1098 -------
1099 mask : `lsst.afw.image.Mask`
1100 Updated mask.
1101 """
1102 if self._doUsePsfMatchedPolygons:
1103 if ("NO_DATA" in altMaskSpans) and ("NO_DATA" in self.config.badMaskPlanes):
1104 # Clear away any other masks outside the validPolygons. These
1105 # pixels are no longer contributing to inexact PSFs, and will
1106 # still be rejected because of NO_DATA.
1107 # self._doUsePsfMatchedPolygons should be True only in
1108 # CompareWarpAssemble. This mask-clearing step must only occur
1109 # *before* applying the new masks below.
1110 for spanSet in altMaskSpans["NO_DATA"]:
1111 spanSet.clippedTo(mask.getBBox()).clearMask(mask, self.getBadPixelMask())
1113 for plane, spanSetList in altMaskSpans.items():
1114 maskClipValue = mask.addMaskPlane(plane)
1115 for spanSet in spanSetList:
1116 spanSet.clippedTo(mask.getBBox()).setMask(mask, 2**maskClipValue)
1117 return mask
1119 def _setValidPolygons(self, warpList, psfMatchedWarpRefList):
1120 """Set the valid polygons for the warps the same as psfMatchedWarps, if
1121 it exists.
1123 Parameters
1124 ----------
1125 warpList : `Iterable` [`lsst.afw.image.Exposure`]
1126 List of warps.
1127 psfMatchedWarpRefList : `Iterable` \
1128 [`lsst.daf.butler.DeferredDatasetHandle`]
1129 List of references to psfMatchedWarps, in the same order as
1130 ``warpList``.
1132 Raises
1133 ------
1134 ValueError
1135 Raised if PSF-matched warps have detectors that are absent in the
1136 (direct) warps.
1137 """
1138 for warp, psfMatchedWarpRef in zip(warpList, psfMatchedWarpRefList):
1139 if psfMatchedWarpRef is None:
1140 continue
1142 psfMatchedCcdTable = psfMatchedWarpRef.get(component="coaddInputs").ccds
1143 ccdTable = warp.getInfo().getCoaddInputs().ccds
1145 # In some (literal) edge cases, a small part of the CCD may be
1146 # present in directWarp that gets excluded in psfMatchedWarp. It is
1147 # okay to leave validPolygon for those CCDs empty. However, the
1148 # converse is not expected, and an error is raised in that case.
1149 if not set(psfMatchedCcdTable["id"]).issubset(ccdTable["id"]):
1150 visit = psfMatchedWarpRef.dataId["visit"]
1151 raise ValueError(f"PSF-matched warp has additional CCDs for {visit=}")
1153 if len(psfMatchedCcdTable) < len(ccdTable):
1154 self.log.debug(
1155 "PSF-matched warp has missing CCDs for visit = %d, which leaves some CCDs in the direct "
1156 "warp without a validPolygon",
1157 psfMatchedWarpRef.dataId["visit"],
1158 )
1160 for psfMatchedCcdRow in psfMatchedCcdTable:
1161 if not psfMatchedCcdRow.validPolygon:
1162 self.log.warning(
1163 "No validPolygon in PSF-matched warp found for %s. This is likely due to a mismatch "
1164 "in the LSST Science Pipelines version used to produce the warps and the current "
1165 "version. To avoid this warning, regenerate the warps with the current version.",
1166 psfMatchedCcdRow.id,
1167 )
1168 else:
1169 ccdRow = ccdTable.find(value=psfMatchedCcdRow.id, key=ccdTable.getIdKey())
1170 ccdRow.validPolygon = psfMatchedCcdRow.validPolygon
1172 def setBrightObjectMasks(self, exposure, brightObjectMasks, dataId=None):
1173 """Set the bright object masks.
1175 Parameters
1176 ----------
1177 exposure : `lsst.afw.image.Exposure`
1178 Exposure under consideration.
1179 brightObjectMasks : `lsst.afw.table`
1180 Table of bright objects to mask.
1181 dataId : `lsst.daf.butler.DataId`, optional
1182 Data identifier dict for patch.
1183 """
1184 if brightObjectMasks is None:
1185 self.log.warning("Unable to apply bright object mask: none supplied")
1186 return
1187 self.log.info("Applying %d bright object masks to %s", len(brightObjectMasks), dataId)
1188 mask = exposure.getMaskedImage().getMask()
1189 wcs = exposure.getWcs()
1190 plateScale = wcs.getPixelScale(exposure.getBBox().getCenter()).asArcseconds()
1192 for rec in brightObjectMasks:
1193 center = geom.PointI(wcs.skyToPixel(rec.getCoord()))
1194 if rec["type"] == "box":
1195 assert rec["angle"] == 0.0, "Angle != 0 for mask object %s" % rec["id"]
1196 width = rec["width"].asArcseconds() / plateScale # convert to pixels
1197 height = rec["height"].asArcseconds() / plateScale # convert to pixels
1199 halfSize = geom.ExtentI(0.5 * width, 0.5 * height)
1200 bbox = geom.Box2I(center - halfSize, center + halfSize)
1202 bbox = geom.BoxI(
1203 geom.PointI(int(center[0] - 0.5 * width), int(center[1] - 0.5 * height)),
1204 geom.PointI(int(center[0] + 0.5 * width), int(center[1] + 0.5 * height)),
1205 )
1206 spans = afwGeom.SpanSet(bbox)
1207 elif rec["type"] == "circle":
1208 radius = int(rec["radius"].asArcseconds() / plateScale) # convert to pixels
1209 spans = afwGeom.SpanSet.fromShape(radius, offset=center)
1210 else:
1211 self.log.warning("Unexpected region type %s at %s", rec["type"], center)
1212 continue
1213 spans.clippedTo(mask.getBBox()).setMask(mask, self.brightObjectBitmask)
1215 def setInexactPsf(self, mask):
1216 """Set INEXACT_PSF mask plane.
1218 If any of the input images isn't represented in the coadd (due to
1219 clipped pixels or chip gaps), the `CoaddPsf` will be inexact. Flag
1220 these pixels.
1222 Parameters
1223 ----------
1224 mask : `lsst.afw.image.Mask`
1225 Coadded exposure's mask, modified in-place.
1226 """
1227 mask.addMaskPlane("INEXACT_PSF")
1228 inexactPsf = mask.getPlaneBitMask("INEXACT_PSF")
1229 sensorEdge = mask.getPlaneBitMask("SENSOR_EDGE") # chip edges (so PSF is discontinuous)
1230 clipped = mask.getPlaneBitMask("CLIPPED") # pixels clipped from coadd
1231 rejected = mask.getPlaneBitMask("REJECTED") # pixels rejected from coadd due to masks
1232 array = mask.getArray()
1233 selected = array & (sensorEdge | clipped | rejected) > 0
1234 array[selected] |= inexactPsf
1236 def filterWarps(self, inputs, goodVisits):
1237 """Return list of only inputRefs with visitId in goodVisits ordered by
1238 goodVisit.
1240 Parameters
1241 ----------
1242 inputs : `list` of `~lsst.pipe.base.connections.DeferredDatasetRef`
1243 List of `lsst.pipe.base.connections.DeferredDatasetRef` with dataId
1244 containing visit.
1245 goodVisit : `dict`
1246 Dictionary with good visitIds as the keys. Value ignored.
1248 Returns
1249 -------
1250 filterInputs : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]
1251 Filtered and sorted list of inputRefs with visitId in goodVisits
1252 ordered by goodVisit.
1253 """
1254 inputWarpDict = {inputRef.ref.dataId["visit"]: inputRef for inputRef in inputs}
1255 filteredInputs = []
1256 for visit in goodVisits.keys():
1257 if visit in inputWarpDict:
1258 filteredInputs.append(inputWarpDict[visit])
1259 return filteredInputs
1262def countMaskFromFootprint(mask, footprint, bitmask, ignoreMask):
1263 """Function to count the number of pixels with a specific mask in a
1264 footprint.
1266 Find the intersection of mask & footprint. Count all pixels in the mask
1267 that are in the intersection that have bitmask set but do not have
1268 ignoreMask set. Return the count.
1270 Parameters
1271 ----------
1272 mask : `lsst.afw.image.Mask`
1273 Mask to define intersection region by.
1274 footprint : `lsst.afw.detection.Footprint`
1275 Footprint to define the intersection region by.
1276 bitmask : `Unknown`
1277 Specific mask that we wish to count the number of occurances of.
1278 ignoreMask : `Unknown`
1279 Pixels to not consider.
1281 Returns
1282 -------
1283 result : `int`
1284 Number of pixels in footprint with specified mask.
1285 """
1286 bbox = footprint.getBBox()
1287 bbox.clip(mask.getBBox(afwImage.PARENT))
1288 fp = afwImage.Mask(bbox)
1289 subMask = mask.Factory(mask, bbox, afwImage.PARENT)
1290 footprint.spans.setMask(fp, bitmask)
1291 return numpy.logical_and(
1292 (subMask.getArray() & fp.getArray()) > 0, (subMask.getArray() & ignoreMask) == 0
1293 ).sum()
1296class CompareWarpAssembleCoaddConnections(AssembleCoaddConnections):
1297 psfMatchedWarps = pipeBase.connectionTypes.Input(
1298 doc=(
1299 "PSF-Matched Warps are required by CompareWarp regardless of the coadd type requested. "
1300 "Only PSF-Matched Warps make sense for image subtraction. "
1301 "Therefore, they must be an additional declared input."
1302 ),
1303 name="{inputCoaddName}Coadd_psfMatchedWarp",
1304 storageClass="ExposureF",
1305 dimensions=("tract", "patch", "skymap", "visit"),
1306 deferLoad=True,
1307 multiple=True,
1308 )
1309 templateCoadd = pipeBase.connectionTypes.Output(
1310 doc=(
1311 "Model of the static sky, used to find temporal artifacts. Typically a PSF-Matched, "
1312 "sigma-clipped coadd. Written if and only if assembleStaticSkyModel.doWrite=True"
1313 ),
1314 name="{outputCoaddName}CoaddPsfMatched",
1315 storageClass="ExposureF",
1316 dimensions=("tract", "patch", "skymap", "band"),
1317 )
1318 artifactMasks = pipeBase.connectionTypes.Output(
1319 doc="Mask of artifacts detected in the coadd",
1320 name="compare_warp_artifact_mask",
1321 storageClass="Mask",
1322 dimensions=("tract", "patch", "skymap", "visit", "instrument"),
1323 multiple=True,
1324 )
1326 def __init__(self, *, config=None):
1327 super().__init__(config=config)
1328 if not config.assembleStaticSkyModel.doWrite:
1329 self.outputs.remove("templateCoadd")
1330 config.validate()
1333class CompareWarpAssembleCoaddConfig(
1334 AssembleCoaddConfig, pipelineConnections=CompareWarpAssembleCoaddConnections
1335):
1336 assembleStaticSkyModel = pexConfig.ConfigurableField(
1337 target=AssembleCoaddTask,
1338 doc="Task to assemble an artifact-free, PSF-matched Coadd to serve as "
1339 "a naive/first-iteration model of the static sky.",
1340 )
1341 detect = pexConfig.ConfigurableField(
1342 target=SourceDetectionTask,
1343 doc="Detect outlier sources on difference between each psfMatched warp and static sky model",
1344 )
1345 detectTemplate = pexConfig.ConfigurableField(
1346 target=SourceDetectionTask,
1347 doc="Detect sources on static sky model. Only used if doPreserveContainedBySource is True",
1348 )
1349 maskStreaks = pexConfig.ConfigurableField(
1350 target=MaskStreaksTask,
1351 doc="Detect streaks on difference between each psfMatched warp and static sky model. Only used if "
1352 "doFilterMorphological is True. Adds a mask plane to an exposure, with the mask plane name set by"
1353 "streakMaskName",
1354 )
1355 streakMaskName = pexConfig.Field(dtype=str, default="STREAK", doc="Name of mask bit used for streaks")
1356 maxNumEpochs = pexConfig.Field(
1357 doc="Charactistic maximum local number of epochs/visits in which an artifact candidate can appear "
1358 "and still be masked. The effective maxNumEpochs is a broken linear function of local "
1359 "number of epochs (N): min(maxFractionEpochsLow*N, maxNumEpochs + maxFractionEpochsHigh*N). "
1360 "For each footprint detected on the image difference between the psfMatched warp and static sky "
1361 "model, if a significant fraction of pixels (defined by spatialThreshold) are residuals in more "
1362 "than the computed effective maxNumEpochs, the artifact candidate is deemed persistant rather "
1363 "than transient and not masked.",
1364 dtype=int,
1365 default=2,
1366 )
1367 maxFractionEpochsLow = pexConfig.RangeField(
1368 doc="Fraction of local number of epochs (N) to use as effective maxNumEpochs for low N. "
1369 "Effective maxNumEpochs = "
1370 "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)",
1371 dtype=float,
1372 default=0.4,
1373 min=0.0,
1374 max=1.0,
1375 )
1376 maxFractionEpochsHigh = pexConfig.RangeField(
1377 doc="Fraction of local number of epochs (N) to use as effective maxNumEpochs for high N. "
1378 "Effective maxNumEpochs = "
1379 "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)",
1380 dtype=float,
1381 default=0.03,
1382 min=0.0,
1383 max=1.0,
1384 )
1385 spatialThreshold = pexConfig.RangeField(
1386 doc="Unitless fraction of pixels defining how much of the outlier region has to meet the "
1387 "temporal criteria. If 0, clip all. If 1, clip none.",
1388 dtype=float,
1389 default=0.5,
1390 min=0.0,
1391 max=1.0,
1392 inclusiveMin=True,
1393 inclusiveMax=True,
1394 )
1395 doScaleWarpVariance = pexConfig.Field(
1396 doc="Rescale Warp variance plane using empirical noise?",
1397 dtype=bool,
1398 default=True,
1399 )
1400 scaleWarpVariance = pexConfig.ConfigurableField(
1401 target=ScaleVarianceTask,
1402 doc="Rescale variance on warps",
1403 )
1404 doPreserveContainedBySource = pexConfig.Field(
1405 doc="Rescue artifacts from clipping that completely lie within a footprint detected"
1406 "on the PsfMatched Template Coadd. Replicates a behavior of SafeClip.",
1407 dtype=bool,
1408 default=True,
1409 )
1410 doPrefilterArtifacts = pexConfig.Field(
1411 doc="Ignore artifact candidates that are mostly covered by the bad pixel mask, "
1412 "because they will be excluded anyway. This prevents them from contributing "
1413 "to the outlier epoch count image and potentially being labeled as persistant."
1414 "'Mostly' is defined by the config 'prefilterArtifactsRatio'.",
1415 dtype=bool,
1416 default=True,
1417 )
1418 prefilterArtifactsMaskPlanes = pexConfig.ListField(
1419 doc="Prefilter artifact candidates that are mostly covered by these bad mask planes.",
1420 dtype=str,
1421 default=("NO_DATA", "BAD", "SUSPECT"),
1422 )
1423 prefilterArtifactsRatio = pexConfig.Field(
1424 doc="Prefilter artifact candidates with less than this fraction overlapping good pixels",
1425 dtype=float,
1426 default=0.05,
1427 )
1428 doFilterMorphological = pexConfig.Field(
1429 doc="Filter artifact candidates based on morphological criteria, i.g. those that appear to "
1430 "be streaks.",
1431 dtype=bool,
1432 default=False,
1433 )
1434 growStreakFp = pexConfig.Field(
1435 doc="Grow streak footprints by this number multiplied by the PSF width", dtype=float, default=5
1436 )
1438 def setDefaults(self):
1439 AssembleCoaddConfig.setDefaults(self)
1440 self.statistic = "MEAN"
1442 # Real EDGE removed by psfMatched NO_DATA border half the width of the
1443 # matching kernel. CompareWarp applies psfMatched EDGE pixels to
1444 # directWarps before assembling.
1445 if "EDGE" in self.badMaskPlanes:
1446 self.badMaskPlanes.remove("EDGE")
1447 self.removeMaskPlanes.append("EDGE")
1448 self.assembleStaticSkyModel.badMaskPlanes = [
1449 "NO_DATA",
1450 ]
1451 self.assembleStaticSkyModel.warpType = "psfMatched"
1452 self.assembleStaticSkyModel.connections.warpType = "psfMatched"
1453 self.assembleStaticSkyModel.statistic = "MEANCLIP"
1454 self.assembleStaticSkyModel.sigmaClip = 2.5
1455 self.assembleStaticSkyModel.clipIter = 3
1456 self.assembleStaticSkyModel.calcErrorFromInputVariance = False
1457 self.assembleStaticSkyModel.doWrite = False
1458 self.detect.doTempLocalBackground = False
1459 self.detect.reEstimateBackground = False
1460 self.detect.returnOriginalFootprints = False
1461 self.detect.thresholdPolarity = "both"
1462 self.detect.thresholdValue = 5
1463 self.detect.minPixels = 4
1464 self.detect.isotropicGrow = True
1465 self.detect.thresholdType = "pixel_stdev"
1466 self.detect.nSigmaToGrow = 0.4
1467 # The default nSigmaToGrow for SourceDetectionTask is already 2.4,
1468 # Explicitly restating because ratio with detect.nSigmaToGrow matters
1469 self.detectTemplate.nSigmaToGrow = 2.4
1470 # Higher thresholds make smaller and fewer protected zones around
1471 # bright stars. Sources with snr < 50 tend to subtract OK empirically
1472 self.detectTemplate.thresholdValue = 50
1473 self.detectTemplate.doTempLocalBackground = False
1474 self.detectTemplate.reEstimateBackground = False
1475 self.detectTemplate.returnOriginalFootprints = False
1477 def validate(self):
1478 super().validate()
1479 if self.assembleStaticSkyModel.doNImage:
1480 raise ValueError(
1481 "No dataset type exists for a PSF-Matched Template N Image."
1482 "Please set assembleStaticSkyModel.doNImage=False"
1483 )
1485 if self.assembleStaticSkyModel.doWrite and (self.warpType == self.assembleStaticSkyModel.warpType):
1486 raise ValueError(
1487 "warpType (%s) == assembleStaticSkyModel.warpType (%s) and will compete for "
1488 "the same dataset name. Please set assembleStaticSkyModel.doWrite to False "
1489 "or warpType to 'direct'. assembleStaticSkyModel.warpType should ways be "
1490 "'PsfMatched'" % (self.warpType, self.assembleStaticSkyModel.warpType)
1491 )
1494class CompareWarpAssembleCoaddTask(AssembleCoaddTask):
1495 """Assemble a compareWarp coadded image from a set of warps
1496 by masking artifacts detected by comparing PSF-matched warps.
1498 In ``AssembleCoaddTask``, we compute the coadd as an clipped mean (i.e.,
1499 we clip outliers). The problem with doing this is that when computing the
1500 coadd PSF at a given location, individual visit PSFs from visits with
1501 outlier pixels contribute to the coadd PSF and cannot be treated correctly.
1502 In this task, we correct for this behavior by creating a new badMaskPlane
1503 'CLIPPED' which marks pixels in the individual warps suspected to contain
1504 an artifact. We populate this plane on the input warps by comparing
1505 PSF-matched warps with a PSF-matched median coadd which serves as a
1506 model of the static sky. Any group of pixels that deviates from the
1507 PSF-matched template coadd by more than config.detect.threshold sigma,
1508 is an artifact candidate. The candidates are then filtered to remove
1509 variable sources and sources that are difficult to subtract such as
1510 bright stars. This filter is configured using the config parameters
1511 ``temporalThreshold`` and ``spatialThreshold``. The temporalThreshold is
1512 the maximum fraction of epochs that the deviation can appear in and still
1513 be considered an artifact. The spatialThreshold is the maximum fraction of
1514 pixels in the footprint of the deviation that appear in other epochs
1515 (where other epochs is defined by the temporalThreshold). If the deviant
1516 region meets this criteria of having a significant percentage of pixels
1517 that deviate in only a few epochs, these pixels have the 'CLIPPED' bit
1518 set in the mask. These regions will not contribute to the final coadd.
1519 Furthermore, any routine to determine the coadd PSF can now be cognizant
1520 of clipped regions. Note that the algorithm implemented by this task is
1521 preliminary and works correctly for HSC data. Parameter modifications and
1522 or considerable redesigning of the algorithm is likley required for other
1523 surveys.
1525 ``CompareWarpAssembleCoaddTask`` sub-classes
1526 ``AssembleCoaddTask`` and instantiates ``AssembleCoaddTask``
1527 as a subtask to generate the TemplateCoadd (the model of the static sky).
1529 Notes
1530 -----
1531 Debugging:
1532 This task supports the following debug variables:
1533 - ``saveCountIm``
1534 If True then save the Epoch Count Image as a fits file in the `figPath`
1535 - ``figPath``
1536 Path to save the debug fits images and figures
1537 """
1539 ConfigClass = CompareWarpAssembleCoaddConfig
1540 _DefaultName = "compareWarpAssembleCoadd"
1542 # See the parent class for docstring.
1543 _doUsePsfMatchedPolygons: bool = True
1545 def __init__(self, *args, **kwargs):
1546 AssembleCoaddTask.__init__(self, *args, **kwargs)
1547 self.makeSubtask("assembleStaticSkyModel")
1548 detectionSchema = afwTable.SourceTable.makeMinimalSchema()
1549 self.makeSubtask("detect", schema=detectionSchema)
1550 if self.config.doPreserveContainedBySource:
1551 self.makeSubtask("detectTemplate", schema=afwTable.SourceTable.makeMinimalSchema())
1552 if self.config.doScaleWarpVariance:
1553 self.makeSubtask("scaleWarpVariance")
1554 if self.config.doFilterMorphological:
1555 self.makeSubtask("maskStreaks")
1557 @utils.inheritDoc(AssembleCoaddTask)
1558 def _makeSupplementaryData(self, butlerQC, inputRefs, outputRefs):
1559 """Generate a templateCoadd to use as a naive model of static sky to
1560 subtract from PSF-Matched warps.
1562 Returns
1563 -------
1564 result : `~lsst.pipe.base.Struct`
1565 Results as a struct with attributes:
1567 ``templateCoadd``
1568 Coadded exposure (`lsst.afw.image.Exposure`).
1569 ``nImage``
1570 Keeps track of exposure count for each pixel
1571 (`lsst.afw.image.ImageU`).
1573 Raises
1574 ------
1575 RuntimeError
1576 Raised if ``templateCoadd`` is `None`.
1577 """
1578 # Ensure that psfMatchedWarps are used as input warps for template
1579 # generation.
1580 staticSkyModelInputRefs = copy.deepcopy(inputRefs)
1581 staticSkyModelInputRefs.inputWarps = inputRefs.psfMatchedWarps
1583 # Because subtasks don't have connections we have to make one.
1584 # The main task's `templateCoadd` is the subtask's `coaddExposure`
1585 staticSkyModelOutputRefs = copy.deepcopy(outputRefs)
1586 if self.config.assembleStaticSkyModel.doWrite:
1587 staticSkyModelOutputRefs.coaddExposure = staticSkyModelOutputRefs.templateCoadd
1588 # Remove template coadd from both subtask's and main tasks outputs,
1589 # because it is handled by the subtask as `coaddExposure`
1590 del outputRefs.templateCoadd
1591 del staticSkyModelOutputRefs.templateCoadd
1593 # A PSF-Matched nImage does not exist as a dataset type
1594 if "nImage" in staticSkyModelOutputRefs.keys():
1595 del staticSkyModelOutputRefs.nImage
1597 templateCoadd = self.assembleStaticSkyModel.runQuantum(
1598 butlerQC, staticSkyModelInputRefs, staticSkyModelOutputRefs
1599 )
1600 if templateCoadd is None:
1601 raise RuntimeError(self._noTemplateMessage(self.assembleStaticSkyModel.warpType))
1603 return pipeBase.Struct(
1604 templateCoadd=templateCoadd.coaddExposure,
1605 nImage=templateCoadd.nImage,
1606 warpRefList=templateCoadd.warpRefList,
1607 imageScalerList=templateCoadd.imageScalerList,
1608 weightList=templateCoadd.weightList,
1609 psfMatchedWarpRefList=inputRefs.psfMatchedWarps,
1610 )
1612 def _noTemplateMessage(self, warpType):
1613 warpName = warpType[0].upper() + warpType[1:]
1614 message = """No %(warpName)s warps were found to build the template coadd which is
1615 required to run CompareWarpAssembleCoaddTask. To continue assembling this type of coadd,
1616 first either rerun makeCoaddTempExp with config.make%(warpName)s=True or
1617 coaddDriver with config.makeCoadTempExp.make%(warpName)s=True, before assembleCoadd.
1619 Alternatively, to use another algorithm with existing warps, retarget the CoaddDriverConfig to
1620 another algorithm like:
1622 from lsst.pipe.tasks.assembleCoadd import SafeClipAssembleCoaddTask
1623 config.assemble.retarget(SafeClipAssembleCoaddTask)
1624 """ % {
1625 "warpName": warpName
1626 }
1627 return message
1629 @utils.inheritDoc(AssembleCoaddTask)
1630 @timeMethod
1631 def run(
1632 self,
1633 skyInfo,
1634 *,
1635 warpRefList,
1636 imageScalerList,
1637 weightList,
1638 psfMatchedWarpRefList,
1639 supplementaryData,
1640 ):
1641 """Notes
1642 -----
1643 Assemble the coadd.
1645 Find artifacts and apply them to the warps' masks creating a list of
1646 alternative masks with a new "CLIPPED" plane and updated "NO_DATA"
1647 plane. Then pass these alternative masks to the base class's ``run``
1648 method.
1649 """
1650 # Check and match the order of the supplementaryData
1651 # (PSF-matched) inputs to the order of the direct inputs,
1652 # so that the artifact mask is applied to the right warp
1654 # For any missing psfMatched warp refs replaced with None,
1655 # findArtifacts() will produce a NO_DATA mask covering the entire bbox.
1656 # As long as NO_DATA is in self.config.badMaskPlanes, these direct
1657 # warps with missing psfMatched warps will not contribute to the coadd.
1658 dataIds = [ref.dataId for ref in warpRefList]
1659 psfMatchedDataIds = [ref.dataId for ref in supplementaryData.warpRefList]
1661 if dataIds != psfMatchedDataIds:
1662 self.log.info("Reordering and or/padding PSF-matched visit input list")
1663 supplementaryData.warpRefList = reorderAndPadList(
1664 supplementaryData.warpRefList, psfMatchedDataIds, dataIds
1665 )
1666 # Remove in DM-49083
1667 supplementaryData.imageScalerList = reorderAndPadList(
1668 supplementaryData.imageScalerList, psfMatchedDataIds, dataIds
1669 )
1671 # Use PSF-Matched Warps (and corresponding scalers) and coadd to find
1672 # artifacts.
1673 spanSetMaskList = self.findArtifacts(
1674 supplementaryData.templateCoadd, supplementaryData.warpRefList, supplementaryData.imageScalerList
1675 )
1677 badMaskPlanes = self.config.badMaskPlanes[:]
1678 badMaskPlanes.append("CLIPPED")
1679 badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes)
1681 result = AssembleCoaddTask.run(
1682 self,
1683 skyInfo,
1684 warpRefList=warpRefList,
1685 imageScalerList=imageScalerList,
1686 weightList=weightList,
1687 altMaskList=spanSetMaskList,
1688 mask=badPixelMask,
1689 psfMatchedWarpRefList=psfMatchedWarpRefList,
1690 )
1692 # Propagate PSF-matched EDGE pixels to coadd SENSOR_EDGE and
1693 # INEXACT_PSF. Psf-Matching moves the real edge inwards.
1694 self.applyAltEdgeMask(result.coaddExposure.maskedImage.mask, spanSetMaskList)
1695 return result
1697 def applyAltEdgeMask(self, mask, altMaskList):
1698 """Propagate alt EDGE mask to SENSOR_EDGE AND INEXACT_PSF planes.
1700 Parameters
1701 ----------
1702 mask : `lsst.afw.image.Mask`
1703 Original mask.
1704 altMaskList : `list` of `dict`
1705 List of Dicts containing ``spanSet`` lists.
1706 Each element contains the new mask plane name (e.g. "CLIPPED
1707 and/or "NO_DATA") as the key, and list of ``SpanSets`` to apply to
1708 the mask.
1709 """
1710 maskValue = mask.getPlaneBitMask(["SENSOR_EDGE", "INEXACT_PSF"])
1711 for visitMask in altMaskList:
1712 if "EDGE" in visitMask:
1713 for spanSet in visitMask["EDGE"]:
1714 spanSet.clippedTo(mask.getBBox()).setMask(mask, maskValue)
1716 def findArtifacts(self, templateCoadd, warpRefList, imageScalerList):
1717 """Find artifacts.
1719 Loop through warps twice. The first loop builds a map with the count
1720 of how many epochs each pixel deviates from the templateCoadd by more
1721 than ``config.chiThreshold`` sigma. The second loop takes each
1722 difference image and filters the artifacts detected in each using
1723 count map to filter out variable sources and sources that are
1724 difficult to subtract cleanly.
1726 Parameters
1727 ----------
1728 templateCoadd : `lsst.afw.image.Exposure`
1729 Exposure to serve as model of static sky.
1730 warpRefList : `list`
1731 List of dataset handles (data references) to warps.
1732 imageScalerList : `list`
1733 List of image scalers.
1734 Deprecated and will be removed after v29.
1736 Returns
1737 -------
1738 altMasks : `list` of `dict`
1739 List of dicts containing information about CLIPPED
1740 (i.e., artifacts), NO_DATA, and EDGE pixels.
1741 """
1742 self.log.debug("Generating Count Image, and mask lists.")
1743 coaddBBox = templateCoadd.getBBox()
1744 slateIm = afwImage.ImageU(coaddBBox)
1745 epochCountImage = afwImage.ImageU(coaddBBox)
1746 nImage = afwImage.ImageU(coaddBBox)
1747 spanSetArtifactList = []
1748 spanSetNoDataMaskList = []
1749 spanSetEdgeList = []
1750 spanSetBadMorphoList = []
1751 badPixelMask = self.getBadPixelMask()
1753 # mask of the warp diffs should = that of only the warp
1754 templateCoadd.mask.clearAllMaskPlanes()
1756 if self.config.doPreserveContainedBySource:
1757 templateFootprints = self.detectTemplate.detectFootprints(templateCoadd)
1758 else:
1759 templateFootprints = None
1761 for warpRef, imageScaler in zip(warpRefList, imageScalerList):
1762 warpDiffExp = self._readAndComputeWarpDiff(warpRef, imageScaler, templateCoadd)
1763 if warpDiffExp is not None:
1764 # This nImage only approximates the final nImage because it
1765 # uses the PSF-matched mask.
1766 nImage.array += (
1767 numpy.isfinite(warpDiffExp.image.array) * ((warpDiffExp.mask.array & badPixelMask) == 0)
1768 ).astype(numpy.uint16)
1769 fpSet = self.detect.detectFootprints(warpDiffExp, doSmooth=False, clearMask=True)
1770 fpSet.positive.merge(fpSet.negative)
1771 footprints = fpSet.positive
1772 slateIm.set(0)
1773 spanSetList = [footprint.spans for footprint in footprints.getFootprints()]
1775 # Remove artifacts due to defects before they contribute to
1776 # the epochCountImage.
1777 if self.config.doPrefilterArtifacts:
1778 spanSetList = self.prefilterArtifacts(spanSetList, warpDiffExp)
1780 # Clear mask before adding prefiltered spanSets
1781 self.detect.clearMask(warpDiffExp.mask)
1782 for spans in spanSetList:
1783 spans.setImage(slateIm, 1, doClip=True)
1784 spans.setMask(warpDiffExp.mask, warpDiffExp.mask.getPlaneBitMask("DETECTED"))
1785 epochCountImage += slateIm
1787 if self.config.doFilterMorphological:
1788 maskName = self.config.streakMaskName
1789 # clear single frame streak mask if it exists already
1790 if maskName in warpDiffExp.mask.getMaskPlaneDict():
1791 warpDiffExp.mask.clearMaskPlane(warpDiffExp.mask.getMaskPlane(maskName))
1792 else:
1793 self.log.debug(f"Did not (need to) clear {maskName} mask because it didn't exist")
1795 _ = self.maskStreaks.run(warpDiffExp)
1796 streakMask = warpDiffExp.mask
1797 initSpanSetStreak = afwGeom.SpanSet.fromMask(
1798 streakMask, streakMask.getPlaneBitMask(maskName)
1799 ).split()
1800 # Pad the streaks to account for low-surface brightness
1801 # wings.
1802 psf = warpDiffExp.getPsf()
1803 spanSetStreak = []
1804 for sset in initSpanSetStreak:
1805 if self.config.doPreserveContainedBySource and templateFootprints is not None:
1806 doKeep = True
1807 for footprint in templateFootprints.positive.getFootprints():
1808 if footprint.spans.contains(sset):
1809 doKeep = False
1810 break
1811 if not doKeep:
1812 continue
1813 psfShape = psf.computeShape(sset.computeCentroid())
1814 dilation = self.config.growStreakFp * psfShape.getDeterminantRadius()
1815 sset_dilated = sset.dilated(int(dilation))
1816 spanSetStreak.append(sset_dilated)
1818 # PSF-Matched warps have less available area (~the matching
1819 # kernel) because the calexps undergo a second convolution.
1820 # Pixels with data in the direct warp but not in the
1821 # PSF-matched warp will not have their artifacts detected.
1822 # NaNs from the PSF-matched warp therefore must be masked in
1823 # the direct warp.
1824 nans = numpy.where(numpy.isnan(warpDiffExp.maskedImage.image.array), 1, 0)
1825 nansMask = afwImage.makeMaskFromArray(nans.astype(afwImage.MaskPixel))
1826 nansMask.setXY0(warpDiffExp.getXY0())
1827 edgeMask = warpDiffExp.mask
1828 spanSetEdgeMask = afwGeom.SpanSet.fromMask(edgeMask, edgeMask.getPlaneBitMask("EDGE")).split()
1829 else:
1830 # If the directWarp has <1% coverage, the psfMatchedWarp can
1831 # have 0% and not exist. In this case, mask the whole epoch.
1832 nansMask = afwImage.MaskX(coaddBBox, 1)
1833 spanSetList = []
1834 spanSetEdgeMask = []
1835 spanSetStreak = []
1837 spanSetNoDataMask = afwGeom.SpanSet.fromMask(nansMask).split()
1839 spanSetNoDataMaskList.append(spanSetNoDataMask)
1840 spanSetArtifactList.append(spanSetList)
1841 spanSetEdgeList.append(spanSetEdgeMask)
1842 if self.config.doFilterMorphological:
1843 spanSetBadMorphoList.append(spanSetStreak)
1845 if lsstDebug.Info(__name__).saveCountIm:
1846 path = self._dataRef2DebugPath("epochCountIm", warpRefList[0], coaddLevel=True)
1847 epochCountImage.writeFits(path)
1849 for i, spanSetList in enumerate(spanSetArtifactList):
1850 if spanSetList:
1851 filteredSpanSetList = self.filterArtifacts(
1852 spanSetList, epochCountImage, nImage, templateFootprints
1853 )
1854 spanSetArtifactList[i] = filteredSpanSetList
1855 if self.config.doFilterMorphological:
1856 spanSetArtifactList[i] += spanSetBadMorphoList[i]
1858 altMasks = []
1859 for artifacts, noData, edge in zip(spanSetArtifactList, spanSetNoDataMaskList, spanSetEdgeList):
1860 altMasks.append({"CLIPPED": artifacts, "NO_DATA": noData, "EDGE": edge})
1861 return altMasks
1863 def prefilterArtifacts(self, spanSetList, exp):
1864 """Remove artifact candidates covered by bad mask plane.
1866 Any future editing of the candidate list that does not depend on
1867 temporal information should go in this method.
1869 Parameters
1870 ----------
1871 spanSetList : `list` [`lsst.afw.geom.SpanSet`]
1872 List of SpanSets representing artifact candidates.
1873 exp : `lsst.afw.image.Exposure`
1874 Exposure containing mask planes used to prefilter.
1876 Returns
1877 -------
1878 returnSpanSetList : `list` [`lsst.afw.geom.SpanSet`]
1879 List of SpanSets with artifacts.
1880 """
1881 badPixelMask = exp.mask.getPlaneBitMask(self.config.prefilterArtifactsMaskPlanes)
1882 goodArr = (exp.mask.array & badPixelMask) == 0
1883 returnSpanSetList = []
1884 bbox = exp.getBBox()
1885 x0, y0 = exp.getXY0()
1886 for i, span in enumerate(spanSetList):
1887 y, x = span.clippedTo(bbox).indices()
1888 yIndexLocal = numpy.array(y) - y0
1889 xIndexLocal = numpy.array(x) - x0
1890 goodRatio = numpy.count_nonzero(goodArr[yIndexLocal, xIndexLocal]) / span.getArea()
1891 if goodRatio > self.config.prefilterArtifactsRatio:
1892 returnSpanSetList.append(span)
1893 return returnSpanSetList
1895 def filterArtifacts(self, spanSetList, epochCountImage, nImage, footprintsToExclude=None):
1896 """Filter artifact candidates.
1898 Parameters
1899 ----------
1900 spanSetList : `list` [`lsst.afw.geom.SpanSet`]
1901 List of SpanSets representing artifact candidates.
1902 epochCountImage : `lsst.afw.image.Image`
1903 Image of accumulated number of warpDiff detections.
1904 nImage : `lsst.afw.image.ImageU`
1905 Image of the accumulated number of total epochs contributing.
1907 Returns
1908 -------
1909 maskSpanSetList : `list` [`lsst.afw.geom.SpanSet`]
1910 List of SpanSets with artifacts.
1911 """
1912 maskSpanSetList = []
1913 x0, y0 = epochCountImage.getXY0()
1914 for i, span in enumerate(spanSetList):
1915 y, x = span.indices()
1916 yIdxLocal = [y1 - y0 for y1 in y]
1917 xIdxLocal = [x1 - x0 for x1 in x]
1918 outlierN = epochCountImage.array[yIdxLocal, xIdxLocal]
1919 totalN = nImage.array[yIdxLocal, xIdxLocal]
1921 # effectiveMaxNumEpochs is broken line (fraction of N) with
1922 # characteristic config.maxNumEpochs.
1923 effMaxNumEpochsHighN = self.config.maxNumEpochs + self.config.maxFractionEpochsHigh * numpy.mean(
1924 totalN
1925 )
1926 effMaxNumEpochsLowN = self.config.maxFractionEpochsLow * numpy.mean(totalN)
1927 effectiveMaxNumEpochs = int(min(effMaxNumEpochsLowN, effMaxNumEpochsHighN))
1928 nPixelsBelowThreshold = numpy.count_nonzero((outlierN > 0) & (outlierN <= effectiveMaxNumEpochs))
1929 percentBelowThreshold = nPixelsBelowThreshold / len(outlierN)
1930 if percentBelowThreshold > self.config.spatialThreshold:
1931 maskSpanSetList.append(span)
1933 if self.config.doPreserveContainedBySource and footprintsToExclude is not None:
1934 # If a candidate is contained by a footprint on the template coadd,
1935 # do not clip.
1936 filteredMaskSpanSetList = []
1937 for span in maskSpanSetList:
1938 doKeep = True
1939 for footprint in footprintsToExclude.positive.getFootprints():
1940 if footprint.spans.contains(span):
1941 doKeep = False
1942 break
1943 if doKeep:
1944 filteredMaskSpanSetList.append(span)
1945 maskSpanSetList = filteredMaskSpanSetList
1947 return maskSpanSetList
1949 def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd):
1950 """Fetch a warp from the butler and return a warpDiff.
1952 Parameters
1953 ----------
1954 warpRef : `lsst.daf.butler.DeferredDatasetHandle`
1955 Dataset handle for the warp.
1956 imageScaler : `lsst.pipe.tasks.scaleZeroPoint.ImageScaler`
1957 An image scaler object.
1958 Deprecated and will be removed after v29 in DM-49083.
1959 templateCoadd : `lsst.afw.image.Exposure`
1960 Exposure to be substracted from the scaled warp.
1962 Returns
1963 -------
1964 warp : `lsst.afw.image.Exposure`
1965 Exposure of the image difference between the warp and template.
1966 """
1967 # If the PSF-Matched warp did not exist for this direct warp
1968 # None is holding its place to maintain order in Gen 3
1969 if warpRef is None:
1970 return None
1972 warp = warpRef.get(parameters={"bbox": templateCoadd.getBBox()})
1973 # direct image scaler OK for PSF-matched Warp.
1974 # Remove in DM-49083.
1975 if imageScaler is not None:
1976 imageScaler.scaleMaskedImage(warp.getMaskedImage())
1977 mi = warp.getMaskedImage()
1978 if self.config.doScaleWarpVariance:
1979 try:
1980 self.scaleWarpVariance.run(mi)
1981 except Exception as exc:
1982 self.log.warning("Unable to rescale variance of warp (%s); leaving it as-is", exc)
1983 mi -= templateCoadd.getMaskedImage()
1984 return warp