Coverage for python / lsst / drp / tasks / assemble_coadd.py: 14%
675 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-31 01:48 -0800
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 = True
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 = 7.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 < 300 tend to subtract OK empirically
1472 self.detectTemplate.thresholdValue = 300
1473 self.detectTemplate.doTempLocalBackground = True
1474 self.detectTemplate.reEstimateBackground = True
1475 self.detectTemplate.background.useApprox = False
1476 self.detectTemplate.returnOriginalFootprints = False
1478 def validate(self):
1479 super().validate()
1480 if self.assembleStaticSkyModel.doNImage:
1481 raise ValueError(
1482 "No dataset type exists for a PSF-Matched Template N Image."
1483 "Please set assembleStaticSkyModel.doNImage=False"
1484 )
1486 if self.assembleStaticSkyModel.doWrite and (self.warpType == self.assembleStaticSkyModel.warpType):
1487 raise ValueError(
1488 "warpType (%s) == assembleStaticSkyModel.warpType (%s) and will compete for "
1489 "the same dataset name. Please set assembleStaticSkyModel.doWrite to False "
1490 "or warpType to 'direct'. assembleStaticSkyModel.warpType should ways be "
1491 "'PsfMatched'" % (self.warpType, self.assembleStaticSkyModel.warpType)
1492 )
1495class CompareWarpAssembleCoaddTask(AssembleCoaddTask):
1496 """Assemble a compareWarp coadded image from a set of warps
1497 by masking artifacts detected by comparing PSF-matched warps.
1499 In ``AssembleCoaddTask``, we compute the coadd as an clipped mean (i.e.,
1500 we clip outliers). The problem with doing this is that when computing the
1501 coadd PSF at a given location, individual visit PSFs from visits with
1502 outlier pixels contribute to the coadd PSF and cannot be treated correctly.
1503 In this task, we correct for this behavior by creating a new badMaskPlane
1504 'CLIPPED' which marks pixels in the individual warps suspected to contain
1505 an artifact. We populate this plane on the input warps by comparing
1506 PSF-matched warps with a PSF-matched median coadd which serves as a
1507 model of the static sky. Any group of pixels that deviates from the
1508 PSF-matched template coadd by more than config.detect.threshold sigma,
1509 is an artifact candidate. The candidates are then filtered to remove
1510 variable sources and sources that are difficult to subtract such as
1511 bright stars. This filter is configured using the config parameters
1512 ``temporalThreshold`` and ``spatialThreshold``. The temporalThreshold is
1513 the maximum fraction of epochs that the deviation can appear in and still
1514 be considered an artifact. The spatialThreshold is the maximum fraction of
1515 pixels in the footprint of the deviation that appear in other epochs
1516 (where other epochs is defined by the temporalThreshold). If the deviant
1517 region meets this criteria of having a significant percentage of pixels
1518 that deviate in only a few epochs, these pixels have the 'CLIPPED' bit
1519 set in the mask. These regions will not contribute to the final coadd.
1520 Furthermore, any routine to determine the coadd PSF can now be cognizant
1521 of clipped regions. Note that the algorithm implemented by this task is
1522 preliminary and works correctly for HSC data. Parameter modifications and
1523 or considerable redesigning of the algorithm is likley required for other
1524 surveys.
1526 ``CompareWarpAssembleCoaddTask`` sub-classes
1527 ``AssembleCoaddTask`` and instantiates ``AssembleCoaddTask``
1528 as a subtask to generate the TemplateCoadd (the model of the static sky).
1530 Notes
1531 -----
1532 Debugging:
1533 This task supports the following debug variables:
1534 - ``saveCountIm``
1535 If True then save the Epoch Count Image as a fits file in the `figPath`
1536 - ``figPath``
1537 Path to save the debug fits images and figures
1538 """
1540 ConfigClass = CompareWarpAssembleCoaddConfig
1541 _DefaultName = "compareWarpAssembleCoadd"
1543 # See the parent class for docstring.
1544 _doUsePsfMatchedPolygons: bool = True
1546 def __init__(self, *args, **kwargs):
1547 AssembleCoaddTask.__init__(self, *args, **kwargs)
1548 self.makeSubtask("assembleStaticSkyModel")
1549 detectionSchema = afwTable.SourceTable.makeMinimalSchema()
1550 self.makeSubtask("detect", schema=detectionSchema)
1551 if self.config.doPreserveContainedBySource:
1552 self.makeSubtask("detectTemplate", schema=afwTable.SourceTable.makeMinimalSchema())
1553 if self.config.doScaleWarpVariance:
1554 self.makeSubtask("scaleWarpVariance")
1555 if self.config.doFilterMorphological:
1556 self.makeSubtask("maskStreaks")
1558 @utils.inheritDoc(AssembleCoaddTask)
1559 def _makeSupplementaryData(self, butlerQC, inputRefs, outputRefs):
1560 """Generate a templateCoadd to use as a naive model of static sky to
1561 subtract from PSF-Matched warps.
1563 Returns
1564 -------
1565 result : `~lsst.pipe.base.Struct`
1566 Results as a struct with attributes:
1568 ``templateCoadd``
1569 Coadded exposure (`lsst.afw.image.Exposure`).
1570 ``nImage``
1571 Keeps track of exposure count for each pixel
1572 (`lsst.afw.image.ImageU`).
1574 Raises
1575 ------
1576 RuntimeError
1577 Raised if ``templateCoadd`` is `None`.
1578 """
1579 # Ensure that psfMatchedWarps are used as input warps for template
1580 # generation.
1581 staticSkyModelInputRefs = copy.deepcopy(inputRefs)
1582 staticSkyModelInputRefs.inputWarps = inputRefs.psfMatchedWarps
1584 # self.assembleStaticSkyModel is not expecting psfMatchedWarps as
1585 # input. But in its runQuantum, ButlerQC would try to be helpful and
1586 # read all of them in to memory.
1587 del staticSkyModelInputRefs.psfMatchedWarps
1589 # Because subtasks don't have connections we have to make one.
1590 # The main task's `templateCoadd` is the subtask's `coaddExposure`
1591 staticSkyModelOutputRefs = copy.deepcopy(outputRefs)
1592 if self.config.assembleStaticSkyModel.doWrite:
1593 staticSkyModelOutputRefs.coaddExposure = staticSkyModelOutputRefs.templateCoadd
1594 # Remove template coadd from both subtask's and main tasks outputs,
1595 # because it is handled by the subtask as `coaddExposure`
1596 del outputRefs.templateCoadd
1597 del staticSkyModelOutputRefs.templateCoadd
1599 # A PSF-Matched nImage does not exist as a dataset type
1600 if "nImage" in staticSkyModelOutputRefs.keys():
1601 del staticSkyModelOutputRefs.nImage
1603 templateCoadd = self.assembleStaticSkyModel.runQuantum(
1604 butlerQC, staticSkyModelInputRefs, staticSkyModelOutputRefs
1605 )
1606 if templateCoadd is None:
1607 raise RuntimeError(self._noTemplateMessage(self.assembleStaticSkyModel.warpType))
1609 return pipeBase.Struct(
1610 templateCoadd=templateCoadd.coaddExposure,
1611 nImage=templateCoadd.nImage,
1612 warpRefList=templateCoadd.warpRefList,
1613 imageScalerList=templateCoadd.imageScalerList,
1614 weightList=templateCoadd.weightList,
1615 psfMatchedWarpRefList=inputRefs.psfMatchedWarps,
1616 )
1618 def _noTemplateMessage(self, warpType):
1619 warpName = warpType[0].upper() + warpType[1:]
1620 message = """No %(warpName)s warps were found to build the template coadd which is
1621 required to run CompareWarpAssembleCoaddTask. To continue assembling this type of coadd,
1622 first either rerun makeCoaddTempExp with config.make%(warpName)s=True or
1623 coaddDriver with config.makeCoadTempExp.make%(warpName)s=True, before assembleCoadd.
1625 Alternatively, to use another algorithm with existing warps, retarget the CoaddDriverConfig to
1626 another algorithm like:
1628 from lsst.pipe.tasks.assembleCoadd import SafeClipAssembleCoaddTask
1629 config.assemble.retarget(SafeClipAssembleCoaddTask)
1630 """ % {
1631 "warpName": warpName
1632 }
1633 return message
1635 @utils.inheritDoc(AssembleCoaddTask)
1636 @timeMethod
1637 def run(
1638 self,
1639 skyInfo,
1640 *,
1641 warpRefList,
1642 imageScalerList,
1643 weightList,
1644 psfMatchedWarpRefList,
1645 supplementaryData,
1646 ):
1647 """Notes
1648 -----
1649 Assemble the coadd.
1651 Find artifacts and apply them to the warps' masks creating a list of
1652 alternative masks with a new "CLIPPED" plane and updated "NO_DATA"
1653 plane. Then pass these alternative masks to the base class's ``run``
1654 method.
1655 """
1656 # Check and match the order of the supplementaryData
1657 # (PSF-matched) inputs to the order of the direct inputs,
1658 # so that the artifact mask is applied to the right warp
1660 # For any missing psfMatched warp refs replaced with None,
1661 # findArtifacts() will produce a NO_DATA mask covering the entire bbox.
1662 # As long as NO_DATA is in self.config.badMaskPlanes, these direct
1663 # warps with missing psfMatched warps will not contribute to the coadd.
1664 dataIds = [ref.dataId for ref in warpRefList]
1665 psfMatchedDataIds = [ref.dataId for ref in supplementaryData.warpRefList]
1667 if dataIds != psfMatchedDataIds:
1668 self.log.info("Reordering and or/padding PSF-matched visit input list")
1669 supplementaryData.warpRefList = reorderAndPadList(
1670 supplementaryData.warpRefList, psfMatchedDataIds, dataIds
1671 )
1672 # Remove in DM-49083
1673 supplementaryData.imageScalerList = reorderAndPadList(
1674 supplementaryData.imageScalerList, psfMatchedDataIds, dataIds
1675 )
1677 # Use PSF-Matched Warps (and corresponding scalers) and coadd to find
1678 # artifacts.
1679 spanSetMaskList = self.findArtifacts(
1680 supplementaryData.templateCoadd, supplementaryData.warpRefList, supplementaryData.imageScalerList
1681 )
1682 # The mask planes are defined globally, so we can load the dict from
1683 # the templateCoadd or the warps.
1684 badMaskPlanes = [
1685 mp
1686 for mp in self.config.badMaskPlanes
1687 if mp in supplementaryData.templateCoadd.mask.getMaskPlaneDict().keys()
1688 ]
1689 badMaskPlanes.append("CLIPPED")
1690 badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes)
1692 result = AssembleCoaddTask.run(
1693 self,
1694 skyInfo,
1695 warpRefList=warpRefList,
1696 imageScalerList=imageScalerList,
1697 weightList=weightList,
1698 altMaskList=spanSetMaskList,
1699 mask=badPixelMask,
1700 psfMatchedWarpRefList=psfMatchedWarpRefList,
1701 )
1703 # Propagate PSF-matched EDGE pixels to coadd SENSOR_EDGE and
1704 # INEXACT_PSF. Psf-Matching moves the real edge inwards.
1705 self.applyAltEdgeMask(result.coaddExposure.maskedImage.mask, spanSetMaskList)
1706 return result
1708 def applyAltEdgeMask(self, mask, altMaskList):
1709 """Propagate alt EDGE mask to SENSOR_EDGE AND INEXACT_PSF planes.
1711 Parameters
1712 ----------
1713 mask : `lsst.afw.image.Mask`
1714 Original mask.
1715 altMaskList : `list` of `dict`
1716 List of Dicts containing ``spanSet`` lists.
1717 Each element contains the new mask plane name (e.g. "CLIPPED
1718 and/or "NO_DATA") as the key, and list of ``SpanSets`` to apply to
1719 the mask.
1720 """
1721 maskValue = mask.getPlaneBitMask(["SENSOR_EDGE", "INEXACT_PSF"])
1722 for visitMask in altMaskList:
1723 if "EDGE" in visitMask:
1724 for spanSet in visitMask["EDGE"]:
1725 spanSet.clippedTo(mask.getBBox()).setMask(mask, maskValue)
1727 def findArtifacts(self, templateCoadd, warpRefList, imageScalerList):
1728 """Find artifacts.
1730 Loop through warps twice. The first loop builds a map with the count
1731 of how many epochs each pixel deviates from the templateCoadd by more
1732 than ``config.chiThreshold`` sigma. The second loop takes each
1733 difference image and filters the artifacts detected in each using
1734 count map to filter out variable sources and sources that are
1735 difficult to subtract cleanly.
1737 Parameters
1738 ----------
1739 templateCoadd : `lsst.afw.image.Exposure`
1740 Exposure to serve as model of static sky.
1741 warpRefList : `list`
1742 List of dataset handles (data references) to warps.
1743 imageScalerList : `list`
1744 List of image scalers.
1745 Deprecated and will be removed after v29.
1747 Returns
1748 -------
1749 altMasks : `list` of `dict`
1750 List of dicts containing information about CLIPPED
1751 (i.e., artifacts), NO_DATA, and EDGE pixels.
1752 """
1753 self.log.debug("Generating Count Image, and mask lists.")
1754 coaddBBox = templateCoadd.getBBox()
1755 slateIm = afwImage.ImageU(coaddBBox)
1756 epochCountImage = afwImage.ImageU(coaddBBox)
1757 nImage = afwImage.ImageU(coaddBBox)
1758 spanSetArtifactList = []
1759 spanSetNoDataMaskList = []
1760 spanSetEdgeList = []
1761 spanSetBadMorphoList = []
1762 badPixelMask = self.getBadPixelMask()
1764 # mask of the warp diffs should = that of only the warp
1765 templateCoadd.mask.clearAllMaskPlanes()
1767 if self.config.doPreserveContainedBySource:
1768 sacrificeToDetection = templateCoadd.clone()
1769 templateFootprints = self.detectTemplate.detectFootprints(sacrificeToDetection)
1770 del sacrificeToDetection
1771 else:
1772 templateFootprints = None
1774 for warpRef, imageScaler in zip(warpRefList, imageScalerList):
1775 warpDiffExp = self._readAndComputeWarpDiff(warpRef, imageScaler, templateCoadd)
1776 if warpDiffExp is not None:
1777 # This nImage only approximates the final nImage because it
1778 # uses the PSF-matched mask.
1779 nImage.array += (
1780 numpy.isfinite(warpDiffExp.image.array) * ((warpDiffExp.mask.array & badPixelMask) == 0)
1781 ).astype(numpy.uint16)
1782 fpSet = self.detect.detectFootprints(warpDiffExp, doSmooth=False, clearMask=True)
1783 fpSet.positive.merge(fpSet.negative)
1784 footprints = fpSet.positive
1785 slateIm.set(0)
1786 spanSetList = [footprint.spans for footprint in footprints.getFootprints()]
1788 # Remove artifacts due to defects before they contribute to
1789 # the epochCountImage.
1790 if self.config.doPrefilterArtifacts:
1791 spanSetList = self.prefilterArtifacts(spanSetList, warpDiffExp)
1793 # Clear mask before adding prefiltered spanSets
1794 self.detect.clearMask(warpDiffExp.mask)
1795 for spans in spanSetList:
1796 spans.setImage(slateIm, 1, doClip=True)
1797 spans.setMask(warpDiffExp.mask, warpDiffExp.mask.getPlaneBitMask("DETECTED"))
1798 epochCountImage += slateIm
1800 if self.config.doFilterMorphological:
1801 maskName = self.config.streakMaskName
1802 # clear single frame streak mask if it exists already
1803 if maskName in warpDiffExp.mask.getMaskPlaneDict():
1804 warpDiffExp.mask.clearMaskPlane(warpDiffExp.mask.getMaskPlane(maskName))
1805 else:
1806 self.log.debug(f"Did not (need to) clear {maskName} mask because it didn't exist")
1808 _ = self.maskStreaks.run(warpDiffExp)
1809 streakMask = warpDiffExp.mask
1810 initSpanSetStreak = afwGeom.SpanSet.fromMask(
1811 streakMask, streakMask.getPlaneBitMask(maskName)
1812 ).split()
1813 # Pad the streaks to account for low-surface brightness
1814 # wings.
1815 psf = warpDiffExp.getPsf()
1816 spanSetStreak = []
1817 for sset in initSpanSetStreak:
1818 if self.config.doPreserveContainedBySource and templateFootprints is not None:
1819 doKeep = True
1820 for footprint in templateFootprints.positive.getFootprints():
1821 if footprint.spans.contains(sset):
1822 doKeep = False
1823 break
1824 if not doKeep:
1825 continue
1826 psfShape = psf.computeShape(sset.computeCentroid())
1827 dilation = self.config.growStreakFp * psfShape.getDeterminantRadius()
1828 sset_dilated = sset.dilated(int(dilation))
1829 spanSetStreak.append(sset_dilated)
1831 # PSF-Matched warps have less available area (~the matching
1832 # kernel) because the calexps undergo a second convolution.
1833 # Pixels with data in the direct warp but not in the
1834 # PSF-matched warp will not have their artifacts detected.
1835 # NaNs from the PSF-matched warp therefore must be masked in
1836 # the direct warp.
1837 nans = numpy.where(numpy.isnan(warpDiffExp.maskedImage.image.array), 1, 0)
1838 nansMask = afwImage.makeMaskFromArray(nans.astype(afwImage.MaskPixel))
1839 nansMask.setXY0(warpDiffExp.getXY0())
1840 edgeMask = warpDiffExp.mask
1841 spanSetEdgeMask = afwGeom.SpanSet.fromMask(edgeMask, edgeMask.getPlaneBitMask("EDGE")).split()
1842 else:
1843 # If the directWarp has <1% coverage, the psfMatchedWarp can
1844 # have 0% and not exist. In this case, mask the whole epoch.
1845 nansMask = afwImage.MaskX(coaddBBox, 1)
1846 spanSetList = []
1847 spanSetEdgeMask = []
1848 spanSetStreak = []
1850 spanSetNoDataMask = afwGeom.SpanSet.fromMask(nansMask).split()
1852 spanSetNoDataMaskList.append(spanSetNoDataMask)
1853 spanSetArtifactList.append(spanSetList)
1854 spanSetEdgeList.append(spanSetEdgeMask)
1855 if self.config.doFilterMorphological:
1856 spanSetBadMorphoList.append(spanSetStreak)
1858 if lsstDebug.Info(__name__).saveCountIm:
1859 path = self._dataRef2DebugPath("epochCountIm", warpRefList[0], coaddLevel=True)
1860 epochCountImage.writeFits(path)
1862 for i, spanSetList in enumerate(spanSetArtifactList):
1863 if spanSetList:
1864 filteredSpanSetList = self.filterArtifacts(
1865 spanSetList, epochCountImage, nImage, templateFootprints
1866 )
1867 spanSetArtifactList[i] = filteredSpanSetList
1868 if self.config.doFilterMorphological:
1869 spanSetArtifactList[i] += spanSetBadMorphoList[i]
1871 altMasks = []
1872 for artifacts, noData, edge in zip(spanSetArtifactList, spanSetNoDataMaskList, spanSetEdgeList):
1873 altMasks.append({"CLIPPED": artifacts, "NO_DATA": noData, "EDGE": edge})
1874 return altMasks
1876 def prefilterArtifacts(self, spanSetList, exp):
1877 """Remove artifact candidates covered by bad mask plane.
1879 Any future editing of the candidate list that does not depend on
1880 temporal information should go in this method.
1882 Parameters
1883 ----------
1884 spanSetList : `list` [`lsst.afw.geom.SpanSet`]
1885 List of SpanSets representing artifact candidates.
1886 exp : `lsst.afw.image.Exposure`
1887 Exposure containing mask planes used to prefilter.
1889 Returns
1890 -------
1891 returnSpanSetList : `list` [`lsst.afw.geom.SpanSet`]
1892 List of SpanSets with artifacts.
1893 """
1894 existingMaskPlanes = [
1895 m for m in self.config.prefilterArtifactsMaskPlanes if m in exp.mask.getMaskPlaneDict()
1896 ]
1897 badPixelMask = exp.mask.getPlaneBitMask(existingMaskPlanes)
1898 goodArr = (exp.mask.array & badPixelMask) == 0
1899 returnSpanSetList = []
1900 bbox = exp.getBBox()
1901 x0, y0 = exp.getXY0()
1902 for i, span in enumerate(spanSetList):
1903 y, x = span.clippedTo(bbox).indices()
1904 yIndexLocal = numpy.array(y) - y0
1905 xIndexLocal = numpy.array(x) - x0
1906 goodRatio = numpy.count_nonzero(goodArr[yIndexLocal, xIndexLocal]) / span.getArea()
1907 if goodRatio > self.config.prefilterArtifactsRatio:
1908 returnSpanSetList.append(span)
1909 return returnSpanSetList
1911 def filterArtifacts(self, spanSetList, epochCountImage, nImage, footprintsToExclude=None):
1912 """Filter artifact candidates.
1914 Parameters
1915 ----------
1916 spanSetList : `list` [`lsst.afw.geom.SpanSet`]
1917 List of SpanSets representing artifact candidates.
1918 epochCountImage : `lsst.afw.image.Image`
1919 Image of accumulated number of warpDiff detections.
1920 nImage : `lsst.afw.image.ImageU`
1921 Image of the accumulated number of total epochs contributing.
1923 Returns
1924 -------
1925 maskSpanSetList : `list` [`lsst.afw.geom.SpanSet`]
1926 List of SpanSets with artifacts.
1927 """
1928 maskSpanSetList = []
1929 x0, y0 = epochCountImage.getXY0()
1930 for i, span in enumerate(spanSetList):
1931 y, x = span.indices()
1932 yIdxLocal = [y1 - y0 for y1 in y]
1933 xIdxLocal = [x1 - x0 for x1 in x]
1934 outlierN = epochCountImage.array[yIdxLocal, xIdxLocal]
1935 totalN = nImage.array[yIdxLocal, xIdxLocal]
1937 # effectiveMaxNumEpochs is broken line (fraction of N) with
1938 # characteristic config.maxNumEpochs.
1939 effMaxNumEpochsHighN = self.config.maxNumEpochs + self.config.maxFractionEpochsHigh * numpy.mean(
1940 totalN
1941 )
1942 effMaxNumEpochsLowN = self.config.maxFractionEpochsLow * numpy.mean(totalN)
1943 effectiveMaxNumEpochs = int(min(effMaxNumEpochsLowN, effMaxNumEpochsHighN))
1944 nPixelsBelowThreshold = numpy.count_nonzero((outlierN > 0) & (outlierN <= effectiveMaxNumEpochs))
1945 percentBelowThreshold = nPixelsBelowThreshold / len(outlierN)
1946 if percentBelowThreshold > self.config.spatialThreshold:
1947 maskSpanSetList.append(span)
1949 if self.config.doPreserveContainedBySource and footprintsToExclude is not None:
1950 # If a candidate is contained by a footprint on the template coadd,
1951 # do not clip.
1952 filteredMaskSpanSetList = []
1953 for span in maskSpanSetList:
1954 doKeep = True
1955 for footprint in footprintsToExclude.positive.getFootprints():
1956 if footprint.spans.contains(span):
1957 doKeep = False
1958 break
1959 if doKeep:
1960 filteredMaskSpanSetList.append(span)
1961 maskSpanSetList = filteredMaskSpanSetList
1963 return maskSpanSetList
1965 def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd):
1966 """Fetch a warp from the butler and return a warpDiff.
1968 Parameters
1969 ----------
1970 warpRef : `lsst.daf.butler.DeferredDatasetHandle`
1971 Dataset handle for the warp.
1972 imageScaler : `lsst.pipe.tasks.scaleZeroPoint.ImageScaler`
1973 An image scaler object.
1974 Deprecated and will be removed after v29 in DM-49083.
1975 templateCoadd : `lsst.afw.image.Exposure`
1976 Exposure to be substracted from the scaled warp.
1978 Returns
1979 -------
1980 warp : `lsst.afw.image.Exposure`
1981 Exposure of the image difference between the warp and template.
1982 """
1983 # If the PSF-Matched warp did not exist for this direct warp
1984 # None is holding its place to maintain order in Gen 3
1985 if warpRef is None:
1986 return None
1988 warp = warpRef.get(parameters={"bbox": templateCoadd.getBBox()})
1989 # direct image scaler OK for PSF-matched Warp.
1990 # Remove in DM-49083.
1991 if imageScaler is not None:
1992 imageScaler.scaleMaskedImage(warp.getMaskedImage())
1993 mi = warp.getMaskedImage()
1994 if self.config.doScaleWarpVariance:
1995 try:
1996 self.scaleWarpVariance.run(mi)
1997 except Exception as exc:
1998 self.log.warning("Unable to rescale variance of warp (%s); leaving it as-is", exc)
1999 mi -= templateCoadd.getMaskedImage()
2000 return warp