Coverage for python / lsst / drp / tasks / assemble_coadd.py: 14%
679 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-21 10:50 +0000
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)
453 try:
454 retStruct = self.run(
455 inputData["skyInfo"],
456 warpRefList=inputs.warpRefList,
457 imageScalerList=inputs.imageScalerList,
458 weightList=inputs.weightList,
459 psfMatchedWarpRefList=inputs.psfMatchedWarpRefList,
460 supplementaryData=supplementaryData,
461 )
462 except pipeBase.AlgorithmError as e:
463 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
464 raise error from e
466 inputData.setdefault("brightObjectMask", None)
467 if self.config.doMaskBrightObjects and inputData["brightObjectMask"] is None:
468 log.warning("doMaskBrightObjects is set to True, but brightObjectMask not loaded")
469 self.processResults(retStruct.coaddExposure, inputData["brightObjectMask"], outputDataId)
471 if self.config.doWriteArtifactMasks:
472 artifactMasksRefList = reorderAndPadList(
473 outputRefs.artifactMasks,
474 [ref.dataId for ref in outputRefs.artifactMasks],
475 [ref.dataId for ref in inputs.warpRefList],
476 )
477 for altMask, warpRef, outputRef in zip(
478 retStruct.altMaskList, inputs.warpRefList, artifactMasksRefList, strict=True
479 ):
480 mask = warpRef.get(component="mask", parameters={"bbox": retStruct.coaddExposure.getBBox()})
481 self.applyAltMaskPlanes(mask, altMask)
482 butlerQC.put(mask, outputRef)
484 if self.config.doWrite:
485 butlerQC.put(retStruct, outputRefs)
487 return retStruct
489 def processResults(self, coaddExposure, brightObjectMasks=None, dataId=None):
490 """Interpolate over missing data and mask bright stars.
492 Parameters
493 ----------
494 coaddExposure : `lsst.afw.image.Exposure`
495 The coadded exposure to process.
496 brightObjectMasks : `lsst.afw.table` or `None`, optional
497 Table of bright objects to mask.
498 dataId : `lsst.daf.butler.DataId` or `None`, optional
499 Data identification.
500 """
501 if self.config.doInterp:
502 self.interpImage.run(coaddExposure.getMaskedImage(), planeName="NO_DATA")
503 # The variance must be positive; work around for DM-3201.
504 varArray = coaddExposure.variance.array
505 with numpy.errstate(invalid="ignore"):
506 varArray[:] = numpy.where(varArray > 0, varArray, numpy.inf)
508 if self.config.doMaskBrightObjects:
509 self.setBrightObjectMasks(coaddExposure, brightObjectMasks, dataId)
511 def _makeSupplementaryData(self, butlerQC, inputRefs, outputRefs):
512 """Make additional inputs to run() specific to subclasses (Gen3).
514 Duplicates interface of `runQuantum` method.
515 Available to be implemented by subclasses only if they need the
516 coadd dataRef/handle for performing preliminary processing before
517 assembling the coadd.
519 Parameters
520 ----------
521 butlerQC : `~lsst.pipe.base.QuantumContext`
522 Gen3 Butler object for fetching additional data products before
523 running the Task specialized for quantum being processed.
524 inputRefs : `~lsst.pipe.base.InputQuantizedConnection`
525 Attributes are the names of the connections describing input
526 dataset types. Values are DatasetRefs that task consumes for the
527 corresponding dataset type. DataIds are guaranteed to match data
528 objects in ``inputData``.
529 outputRefs : `~lsst.pipe.base.OutputQuantizedConnection`
530 Attributes are the names of the connections describing output
531 dataset types. Values are DatasetRefs that task is to produce
532 for the corresponding dataset type.
533 """
534 return pipeBase.Struct()
536 def prepareInputs(self, refList, coadd_bbox, psfMatchedWarpRefList=None):
537 """Prepare the input warps for coaddition by measuring the weight for
538 each warp.
540 Before coadding these Warps together compute the weight for each
541 Warp.
543 Parameters
544 ----------
545 refList : `list`
546 List of dataset handles (data references) to warp.
547 psfMatchedWarpRefList : `list` | None, optional
548 List of dataset handles (data references) to psfMatchedWarp.
550 Returns
551 -------
552 result : `~lsst.pipe.base.Struct`
553 Results as a struct with attributes:
555 ``warpRefList``
556 `list` of dataset handles (data references) to warp.
557 ``weightList``
558 `list` of weightings.
559 ``imageScalerList``
560 `list` of image scalers.
561 Deprecated and will be removed in DM-49083.
562 """
563 statsCtrl = afwMath.StatisticsControl()
564 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
565 statsCtrl.setNumIter(self.config.clipIter)
566 statsCtrl.setAndMask(self.getBadPixelMask())
567 statsCtrl.setNanSafe(True)
568 # compute warpRefList: a list of warpRef that actually exist
569 # and weightList: a list of the weight of the associated coadd warp
570 # and imageScalerList: a list of scale factors for the associated coadd
571 # warp. (deprecated and will be removed in DM-49083)
572 warpRefList = []
573 weightList = []
574 # Remove in DM-49083
575 imageScalerList = []
576 outputPsfMatchedWarpRefList = []
578 # Convert psfMatchedWarpRefList to a dict, so we can look them up by
579 # their dataId.
580 # Note: Some warps may not have a corresponding psfMatchedWarp, which
581 # could have happened due to a failure in the PSF matching, or not
582 # having enough good pixels to support the PSF-matching kernel.
584 if psfMatchedWarpRefList is None:
585 psfMatchedWarpRefDict = {ref.dataId: None for ref in refList}
586 else:
587 psfMatchedWarpRefDict = {ref.dataId: ref for ref in psfMatchedWarpRefList}
589 warpName = self.getTempExpDatasetName(self.warpType)
590 for warpRef in refList:
591 warp = warpRef.get(parameters={"bbox": coadd_bbox})
592 # Ignore any input warp that is empty of data
593 if numpy.isnan(warp.image.array).all():
594 continue
595 maskedImage = warp.getMaskedImage()
597 # Allow users to scale the zero point for backward
598 # compatibility. Remove imageScalarList in DM-49083.
599 if self.config.doScaleZeroPoint:
600 imageScaler = self.scaleZeroPoint.computeImageScaler(
601 exposure=warp,
602 dataRef=warpRef, # FIXME
603 )
604 try:
605 imageScaler.scaleMaskedImage(maskedImage)
606 except Exception as e:
607 self.log.warning("Scaling failed for %s (skipping it): %s", warpRef.dataId, e)
608 continue
609 imageScalerList.append(imageScaler)
610 else:
611 imageScalerList.append(None)
612 if "BUNIT" not in warp.metadata:
613 raise ValueError(f"Warp {warpRef.dataId} has no BUNIT metadata")
614 if warp.metadata["BUNIT"] != "nJy":
615 raise ValueError(
616 f"Warp {warpRef.dataId} has BUNIT {warp.metadata['BUNIT']}, expected nJy"
617 )
618 statObj = afwMath.makeStatistics(
619 maskedImage.getVariance(), maskedImage.getMask(), afwMath.MEANCLIP, statsCtrl
620 )
621 meanVar, meanVarErr = statObj.getResult(afwMath.MEANCLIP)
622 if meanVar > 0.0 and numpy.isfinite(meanVar):
623 weight = 1.0 / float(meanVar)
624 self.log.info("Weight of %s %s = %0.3f", warpName, warpRef.dataId, weight)
625 else:
626 self.log.warning("Non-finite weight for %s: skipping", warpRef.dataId)
627 continue
629 del maskedImage
630 del warp
632 warpRefList.append(warpRef)
633 weightList.append(weight)
635 dataId = warpRef.dataId
636 psfMatchedWarpRef = psfMatchedWarpRefDict.get(dataId, None)
637 outputPsfMatchedWarpRefList.append(psfMatchedWarpRef)
639 return pipeBase.Struct(
640 warpRefList=warpRefList,
641 weightList=weightList,
642 imageScalerList=imageScalerList,
643 psfMatchedWarpRefList=outputPsfMatchedWarpRefList,
644 )
646 def prepareStats(self, mask=None):
647 """Prepare the statistics for coadding images.
649 Parameters
650 ----------
651 mask : `int`, optional
652 Bit mask value to exclude from coaddition.
654 Returns
655 -------
656 stats : `~lsst.pipe.base.Struct`
657 Statistics as a struct with attributes:
659 ``statsCtrl``
660 Statistics control object for coadd
661 (`~lsst.afw.math.StatisticsControl`).
662 ``statsFlags``
663 Statistic for coadd (`~lsst.afw.math.Property`).
664 """
665 if mask is None:
666 mask = self.getBadPixelMask()
667 statsCtrl = afwMath.StatisticsControl()
668 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
669 statsCtrl.setNumIter(self.config.clipIter)
670 statsCtrl.setAndMask(mask)
671 statsCtrl.setNanSafe(True)
672 statsCtrl.setWeighted(True)
673 statsCtrl.setCalcErrorFromInputVariance(self.config.calcErrorFromInputVariance)
674 for plane, threshold in self.config.maskPropagationThresholds.items():
675 bit = afwImage.Mask.getMaskPlane(plane)
676 statsCtrl.setMaskPropagationThreshold(bit, threshold)
677 statsFlags = afwMath.stringToStatisticsProperty(self.config.statistic)
678 return pipeBase.Struct(ctrl=statsCtrl, flags=statsFlags)
680 @timeMethod
681 def run(
682 self,
683 skyInfo,
684 *,
685 warpRefList,
686 imageScalerList,
687 weightList,
688 psfMatchedWarpRefList=None,
689 altMaskList=None,
690 mask=None,
691 supplementaryData=None,
692 ):
693 """Assemble a coadd from input warps.
695 Assemble the coadd using the provided list of coaddTempExps. Since
696 the full coadd covers a patch (a large area), the assembly is
697 performed over small areas on the image at a time in order to
698 conserve memory usage. Iterate over subregions within the outer
699 bbox of the patch using `assembleSubregion` to stack the corresponding
700 subregions from the coaddTempExps with the statistic specified.
701 Set the edge bits the coadd mask based on the weight map.
703 Parameters
704 ----------
705 skyInfo : `~lsst.pipe.base.Struct`
706 Struct with geometric information about the patch.
707 warpRefList : `list`
708 List of dataset handles (data references) to Warps
709 (previously called CoaddTempExps).
710 imageScalerList : `list`
711 List of image scalers. Deprecated and will be removed after v29
712 in DM-49083.
713 weightList : `list`
714 List of weights.
715 psfMatchedWarpRefList : `list`, optional
716 List of dataset handles (data references) to psfMatchedWarps.
717 altMaskList : `list`, optional
718 List of alternate masks to use rather than those stored with
719 warp.
720 mask : `int`, optional
721 Bit mask value to exclude from coaddition.
722 supplementaryData : `~lsst.pipe.base.Struct`, optional
723 Struct with additional data products needed to assemble coadd.
724 Only used by subclasses that implement ``_makeSupplementaryData``
725 and override `run`.
727 Returns
728 -------
729 result : `~lsst.pipe.base.Struct`
730 Results as a struct with attributes:
732 ``coaddExposure``
733 Coadded exposure (`~lsst.afw.image.Exposure`).
734 ``nImage``
735 Exposure count image (`~lsst.afw.image.Image`), if requested.
736 ``inputMap``
737 Bit-wise map of inputs, if requested.
738 ``warpRefList``
739 Input list of dataset handles (data refs) to the warps
740 (`~lsst.daf.butler.DeferredDatasetHandle`) (unmodified).
741 ``imageScalerList``
742 Input list of image scalers (`list`) (unmodified).
743 Deprecated and will be removed after v29 in DM-49083.
744 ``weightList``
745 Input list of weights (`list`) (unmodified).
747 Raises
748 ------
749 lsst.pipe.base.NoWorkFound
750 Raised if no dataset handles (data references) are provided.
751 """
752 warpName = self.getTempExpDatasetName(self.warpType)
753 self.log.info("Assembling %s %s", len(warpRefList), warpName)
754 if not warpRefList:
755 raise pipeBase.NoWorkFound("No exposures provided for co-addition.")
757 stats = self.prepareStats(mask=mask)
759 if altMaskList is None:
760 altMaskList = [None] * len(warpRefList)
762 coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs)
763 # Deprecated, keep only the `else` branch in DM-49083.
764 if self.config.doScaleZeroPoint:
765 coaddExposure.setPhotoCalib(self.scaleZeroPoint.getPhotoCalib())
766 else:
767 coaddExposure.setPhotoCalib(afwImage.PhotoCalib(1.0))
768 coaddExposure.getInfo().setCoaddInputs(self.inputRecorder.makeCoaddInputs())
769 self.assembleMetadata(coaddExposure, warpRefList, weightList, psfMatchedWarpRefList)
770 coaddMaskedImage = coaddExposure.getMaskedImage()
771 subregionSizeArr = self.config.subregionSize
772 subregionSize = geom.Extent2I(subregionSizeArr[0], subregionSizeArr[1])
773 # if nImage is requested, create a zero one which can be passed to
774 # assembleSubregion.
775 if self.config.doNImage:
776 nImage = afwImage.ImageU(skyInfo.bbox)
777 else:
778 nImage = None
779 # If inputMap is requested, create the initial version that can be
780 # masked in assembleSubregion.
781 if self.config.doInputMap:
782 self.inputMapper.build_ccd_input_map(
783 skyInfo.bbox, skyInfo.wcs, coaddExposure.getInfo().getCoaddInputs().ccds
784 )
786 if self.config.doOnlineForMean and self.config.statistic == "MEAN":
787 try:
788 self.assembleOnlineMeanCoadd(
789 coaddExposure,
790 warpRefList,
791 imageScalerList,
792 weightList,
793 altMaskList,
794 stats.ctrl,
795 nImage=nImage,
796 )
797 except Exception as e:
798 self.log.exception("Cannot compute online coadd %s", e)
799 raise
800 else:
801 for subBBox in subBBoxIter(skyInfo.bbox, subregionSize):
802 try:
803 self.assembleSubregion(
804 coaddExposure,
805 subBBox,
806 warpRefList,
807 imageScalerList,
808 weightList,
809 altMaskList,
810 stats.flags,
811 stats.ctrl,
812 nImage=nImage,
813 )
814 except Exception as e:
815 self.log.exception("Cannot compute coadd %s: %s", subBBox, e)
816 raise
818 # If inputMap is requested, we must finalize the map after the
819 # accumulation.
820 if self.config.doInputMap:
821 self.inputMapper.finalize_ccd_input_map_mask()
822 inputMap = self.inputMapper.ccd_input_map
823 else:
824 inputMap = None
826 self.setInexactPsf(coaddMaskedImage.getMask())
827 # Despite the name, the following doesn't really deal with "EDGE"
828 # pixels: it identifies pixels that didn't receive any unmasked inputs
829 # (as occurs around the edge of the field).
830 coaddUtils.setCoaddEdgeBits(coaddMaskedImage.getMask(), coaddMaskedImage.getVariance())
831 return pipeBase.Struct(
832 coaddExposure=coaddExposure,
833 nImage=nImage,
834 warpRefList=warpRefList,
835 imageScalerList=imageScalerList,
836 weightList=weightList,
837 inputMap=inputMap,
838 altMaskList=altMaskList,
839 )
841 def assembleMetadata(self, coaddExposure, warpRefList, weightList, psfMatchedWarpRefList=None):
842 """Set the metadata for the coadd.
844 This basic implementation sets the filter from the first input.
846 Parameters
847 ----------
848 coaddExposure : `lsst.afw.image.Exposure`
849 The target exposure for the coadd.
850 warpRefList : `list`
851 List of dataset handles (data references) to warp.
852 weightList : `list`
853 List of weights.
854 psfMatchedWarpRefList : `list` | None, optional
855 List of dataset handles (data references) to psfMatchedWarps.
857 Raises
858 ------
859 AssertionError
860 Raised if there is a length mismatch.
861 """
862 assert len(warpRefList) == len(weightList), "Length mismatch"
864 if psfMatchedWarpRefList:
865 assert len(warpRefList) == len(psfMatchedWarpRefList), "Length mismatch"
867 # We load a single pixel of each coaddTempExp, because we just want to
868 # get at the metadata (and we need more than just the PropertySet that
869 # contains the header), which is not possible with the current butler
870 # (see #2777).
871 bbox = geom.Box2I(coaddExposure.getBBox().getMin(), geom.Extent2I(1, 1))
873 warpList = [warpRef.get(parameters={"bbox": bbox}) for warpRef in warpRefList]
875 numCcds = sum(len(warp.getInfo().getCoaddInputs().ccds) for warp in warpList)
877 # Set the coadd FilterLabel to the band of the first input exposure:
878 # Coadds are calibrated, so the physical label is now meaningless.
879 coaddExposure.setFilter(afwImage.FilterLabel(warpList[0].getFilter().bandLabel))
880 coaddInputs = coaddExposure.getInfo().getCoaddInputs()
881 coaddInputs.ccds.reserve(numCcds)
882 coaddInputs.visits.reserve(len(warpList))
884 # Set the exposure units to nJy
885 # No need to check after DM-49083.
886 if not self.config.doScaleZeroPoint:
887 coaddExposure.metadata["BUNIT"] = "nJy"
889 # psfMatchedWarpRefList should be empty except in CompareWarpCoadd.
890 if self._doUsePsfMatchedPolygons:
891 # Set validPolygons for warp before addVisitToCoadd
892 self._setValidPolygons(warpList, psfMatchedWarpRefList)
894 for warp, weight in zip(warpList, weightList):
895 self.inputRecorder.addVisitToCoadd(coaddInputs, warp, weight)
897 coaddInputs.visits.sort()
898 coaddInputs.ccds.sort()
899 if self.warpType == "psfMatched":
900 # The modelPsf BBox for a psfMatchedWarp/coaddTempExp was
901 # dynamically defined by ModelPsfMatchTask as the square box
902 # bounding its spatially-variable, pre-matched WarpedPsf.
903 # Likewise, set the PSF of a PSF-Matched Coadd to the modelPsf
904 # having the maximum width (sufficient because square)
905 modelPsfList = [warp.getPsf() for warp in warpList]
906 modelPsfWidthList = [
907 modelPsf.computeBBox(modelPsf.getAveragePosition()).getWidth() for modelPsf in modelPsfList
908 ]
909 psf = modelPsfList[modelPsfWidthList.index(max(modelPsfWidthList))]
910 else:
911 psf = measAlg.CoaddPsf(
912 coaddInputs.ccds, coaddExposure.getWcs(), self.config.coaddPsf.makeControl()
913 )
914 coaddExposure.setPsf(psf)
915 apCorrMap = measAlg.makeCoaddApCorrMap(
916 coaddInputs.ccds, coaddExposure.getBBox(afwImage.PARENT), coaddExposure.getWcs()
917 )
918 coaddExposure.getInfo().setApCorrMap(apCorrMap)
919 if self.config.doAttachTransmissionCurve:
920 transmissionCurve = measAlg.makeCoaddTransmissionCurve(coaddExposure.getWcs(), coaddInputs.ccds)
921 coaddExposure.getInfo().setTransmissionCurve(transmissionCurve)
923 def assembleSubregion(
924 self,
925 coaddExposure,
926 bbox,
927 warpRefList,
928 imageScalerList,
929 weightList,
930 altMaskList,
931 statsFlags,
932 statsCtrl,
933 nImage=None,
934 ):
935 """Assemble the coadd for a sub-region.
937 For each coaddTempExp, check for (and swap in) an alternative mask
938 if one is passed. Remove mask planes listed in
939 `config.removeMaskPlanes`. Finally, stack the actual exposures using
940 `lsst.afw.math.statisticsStack` with the statistic specified by
941 statsFlags. Typically, the statsFlag will be one of lsst.afw.math.MEAN
942 for a mean-stack or `lsst.afw.math.MEANCLIP` for outlier rejection
943 using an N-sigma clipped mean where N and iterations are specified by
944 statsCtrl. Assign the stacked subregion back to the coadd.
946 Parameters
947 ----------
948 coaddExposure : `lsst.afw.image.Exposure`
949 The target exposure for the coadd.
950 bbox : `lsst.geom.Box`
951 Sub-region to coadd.
952 warpRefList : `list`
953 List of dataset handles (data references) to warp.
954 imageScalerList : `list`
955 List of image scalers.
956 Deprecated and will be removed after v29 in DM-49083.
957 weightList : `list`
958 List of weights.
959 altMaskList : `list`
960 List of alternate masks to use rather than those stored with
961 warp, or None. Each element is dict with keys = mask plane
962 name to which to add the spans.
963 statsFlags : `lsst.afw.math.Property`
964 Property object for statistic for coadd.
965 statsCtrl : `lsst.afw.math.StatisticsControl`
966 Statistics control object for coadd.
967 nImage : `lsst.afw.image.ImageU`, optional
968 Keeps track of exposure count for each pixel.
969 """
970 self.log.info("Stacking subregion %s", bbox)
972 coaddExposure.mask.addMaskPlane("REJECTED")
973 coaddExposure.mask.addMaskPlane("CLIPPED")
974 coaddExposure.mask.addMaskPlane("SENSOR_EDGE")
975 maskMap = setRejectedMaskMapping(statsCtrl)
976 clipped = afwImage.Mask.getPlaneBitMask("CLIPPED")
977 maskedImageList = []
978 if nImage is not None:
979 subNImage = afwImage.ImageU(bbox.getWidth(), bbox.getHeight())
980 for warpRef, imageScaler, altMask in zip(warpRefList, imageScalerList, altMaskList):
981 exposure = warpRef.get(parameters={"bbox": bbox})
983 maskedImage = exposure.getMaskedImage()
984 mask = maskedImage.getMask()
985 if altMask is not None:
986 self.applyAltMaskPlanes(mask, altMask)
987 # Remove in DM-49083
988 if imageScaler is not None:
989 imageScaler.scaleMaskedImage(maskedImage)
991 # Add 1 for each pixel which is not excluded by the exclude mask.
992 # In legacyCoadd, pixels may also be excluded by
993 # afwMath.statisticsStack.
994 if nImage is not None:
995 subNImage.getArray()[maskedImage.getMask().getArray() & statsCtrl.getAndMask() == 0] += 1
996 if self.config.removeMaskPlanes:
997 removeMaskPlanes(maskedImage.mask, self.config.removeMaskPlanes, logger=self.log)
998 maskedImageList.append(maskedImage)
1000 if self.config.doInputMap:
1001 visit = exposure.getInfo().getCoaddInputs().visits[0].getId()
1002 self.inputMapper.mask_warp_bbox(bbox, visit, mask, statsCtrl.getAndMask())
1004 with self.timer("stack"):
1005 coaddSubregion = afwMath.statisticsStack(
1006 maskedImageList,
1007 statsFlags,
1008 statsCtrl,
1009 weightList,
1010 clipped, # also set output to CLIPPED if sigma-clipped
1011 maskMap,
1012 )
1013 coaddExposure.maskedImage.assign(coaddSubregion, bbox)
1014 if nImage is not None:
1015 nImage.assign(subNImage, bbox)
1017 def assembleOnlineMeanCoadd(
1018 self, coaddExposure, warpRefList, imageScalerList, weightList, altMaskList, statsCtrl, nImage=None
1019 ):
1020 """Assemble the coadd using the "online" method.
1022 This method takes a running sum of images and weights to save memory.
1023 It only works for MEAN statistics.
1025 Parameters
1026 ----------
1027 coaddExposure : `lsst.afw.image.Exposure`
1028 The target exposure for the coadd.
1029 warpRefList : `list`
1030 List of dataset handles (data references) to warp.
1031 imageScalerList : `list`
1032 List of image scalers.
1033 Deprecated and will be removed after v29 in DM-49083.
1034 weightList : `list`
1035 List of weights.
1036 altMaskList : `list`
1037 List of alternate masks to use rather than those stored with
1038 warp, or None. Each element is dict with keys = mask plane
1039 name to which to add the spans.
1040 statsCtrl : `lsst.afw.math.StatisticsControl`
1041 Statistics control object for coadd.
1042 nImage : `lsst.afw.image.ImageU`, optional
1043 Keeps track of exposure count for each pixel.
1044 """
1045 self.log.debug("Computing online coadd.")
1047 coaddExposure.mask.addMaskPlane("REJECTED")
1048 coaddExposure.mask.addMaskPlane("CLIPPED")
1049 coaddExposure.mask.addMaskPlane("SENSOR_EDGE")
1050 maskMap = setRejectedMaskMapping(statsCtrl)
1051 thresholdDict = AccumulatorMeanStack.stats_ctrl_to_threshold_dict(statsCtrl)
1053 bbox = coaddExposure.maskedImage.getBBox()
1055 stacker = AccumulatorMeanStack(
1056 coaddExposure.image.array.shape,
1057 statsCtrl.getAndMask(),
1058 mask_threshold_dict=thresholdDict,
1059 mask_map=maskMap,
1060 no_good_pixels_mask=statsCtrl.getNoGoodPixelsMask(),
1061 calc_error_from_input_variance=self.config.calcErrorFromInputVariance,
1062 compute_n_image=(nImage is not None),
1063 )
1065 for warpRef, imageScaler, altMask, weight in zip(
1066 warpRefList, imageScalerList, altMaskList, weightList
1067 ):
1068 exposure = warpRef.get(parameters={"bbox": bbox})
1069 maskedImage = exposure.getMaskedImage()
1070 mask = maskedImage.getMask()
1071 if altMask is not None:
1072 self.applyAltMaskPlanes(mask, altMask)
1073 # Remove in DM-49083.
1074 if imageScaler is not None:
1075 imageScaler.scaleMaskedImage(maskedImage)
1076 if self.config.removeMaskPlanes:
1077 removeMaskPlanes(maskedImage.mask, self.config.removeMaskPlanes, logger=self.log)
1079 stacker.add_masked_image(maskedImage, weight=weight)
1081 if self.config.doInputMap:
1082 visit = exposure.getInfo().getCoaddInputs().visits[0].getId()
1083 self.inputMapper.mask_warp_bbox(bbox, visit, mask, statsCtrl.getAndMask())
1085 stacker.fill_stacked_masked_image(coaddExposure.maskedImage)
1087 if nImage is not None:
1088 nImage.array[:, :] = stacker.n_image
1090 def applyAltMaskPlanes(self, mask, altMaskSpans):
1091 """Apply in place alt mask formatted as SpanSets to a mask.
1093 Parameters
1094 ----------
1095 mask : `lsst.afw.image.Mask`
1096 Original mask.
1097 altMaskSpans : `dict`
1098 SpanSet lists to apply. Each element contains the new mask
1099 plane name (e.g. "CLIPPED and/or "NO_DATA") as the key,
1100 and list of SpanSets to apply to the mask.
1102 Returns
1103 -------
1104 mask : `lsst.afw.image.Mask`
1105 Updated mask.
1106 """
1107 if self._doUsePsfMatchedPolygons:
1108 if ("NO_DATA" in altMaskSpans) and ("NO_DATA" in self.config.badMaskPlanes):
1109 # Clear away any other masks outside the validPolygons. These
1110 # pixels are no longer contributing to inexact PSFs, and will
1111 # still be rejected because of NO_DATA.
1112 # self._doUsePsfMatchedPolygons should be True only in
1113 # CompareWarpAssemble. This mask-clearing step must only occur
1114 # *before* applying the new masks below.
1115 for spanSet in altMaskSpans["NO_DATA"]:
1116 spanSet.clippedTo(mask.getBBox()).clearMask(mask, self.getBadPixelMask())
1118 for plane, spanSetList in altMaskSpans.items():
1119 maskClipValue = mask.addMaskPlane(plane)
1120 for spanSet in spanSetList:
1121 spanSet.clippedTo(mask.getBBox()).setMask(mask, 2**maskClipValue)
1122 return mask
1124 def _setValidPolygons(self, warpList, psfMatchedWarpRefList):
1125 """Set the valid polygons for the warps the same as psfMatchedWarps, if
1126 it exists.
1128 Parameters
1129 ----------
1130 warpList : `Iterable` [`lsst.afw.image.Exposure`]
1131 List of warps.
1132 psfMatchedWarpRefList : `Iterable` \
1133 [`lsst.daf.butler.DeferredDatasetHandle`]
1134 List of references to psfMatchedWarps, in the same order as
1135 ``warpList``.
1137 Raises
1138 ------
1139 ValueError
1140 Raised if PSF-matched warps have detectors that are absent in the
1141 (direct) warps.
1142 """
1143 for warp, psfMatchedWarpRef in zip(warpList, psfMatchedWarpRefList):
1144 if psfMatchedWarpRef is None:
1145 continue
1147 psfMatchedCcdTable = psfMatchedWarpRef.get(component="coaddInputs").ccds
1148 ccdTable = warp.getInfo().getCoaddInputs().ccds
1150 # In some (literal) edge cases, a small part of the CCD may be
1151 # present in directWarp that gets excluded in psfMatchedWarp. It is
1152 # okay to leave validPolygon for those CCDs empty. However, the
1153 # converse is not expected, and an error is raised in that case.
1154 if not set(psfMatchedCcdTable["id"]).issubset(ccdTable["id"]):
1155 visit = psfMatchedWarpRef.dataId["visit"]
1156 raise ValueError(f"PSF-matched warp has additional CCDs for {visit=}")
1158 if len(psfMatchedCcdTable) < len(ccdTable):
1159 self.log.debug(
1160 "PSF-matched warp has missing CCDs for visit = %d, which leaves some CCDs in the direct "
1161 "warp without a validPolygon",
1162 psfMatchedWarpRef.dataId["visit"],
1163 )
1165 for psfMatchedCcdRow in psfMatchedCcdTable:
1166 if not psfMatchedCcdRow.validPolygon:
1167 self.log.warning(
1168 "No validPolygon in PSF-matched warp found for %s. This is likely due to a mismatch "
1169 "in the LSST Science Pipelines version used to produce the warps and the current "
1170 "version. To avoid this warning, regenerate the warps with the current version.",
1171 psfMatchedCcdRow.id,
1172 )
1173 else:
1174 ccdRow = ccdTable.find(value=psfMatchedCcdRow.id, key=ccdTable.getIdKey())
1175 ccdRow.validPolygon = psfMatchedCcdRow.validPolygon
1177 def setBrightObjectMasks(self, exposure, brightObjectMasks, dataId=None):
1178 """Set the bright object masks.
1180 Parameters
1181 ----------
1182 exposure : `lsst.afw.image.Exposure`
1183 Exposure under consideration.
1184 brightObjectMasks : `lsst.afw.table`
1185 Table of bright objects to mask.
1186 dataId : `lsst.daf.butler.DataId`, optional
1187 Data identifier dict for patch.
1188 """
1189 if brightObjectMasks is None:
1190 self.log.warning("Unable to apply bright object mask: none supplied")
1191 return
1192 self.log.info("Applying %d bright object masks to %s", len(brightObjectMasks), dataId)
1193 mask = exposure.getMaskedImage().getMask()
1194 wcs = exposure.getWcs()
1195 plateScale = wcs.getPixelScale(exposure.getBBox().getCenter()).asArcseconds()
1197 for rec in brightObjectMasks:
1198 center = geom.PointI(wcs.skyToPixel(rec.getCoord()))
1199 if rec["type"] == "box":
1200 assert rec["angle"] == 0.0, "Angle != 0 for mask object %s" % rec["id"]
1201 width = rec["width"].asArcseconds() / plateScale # convert to pixels
1202 height = rec["height"].asArcseconds() / plateScale # convert to pixels
1204 halfSize = geom.ExtentI(0.5 * width, 0.5 * height)
1205 bbox = geom.Box2I(center - halfSize, center + halfSize)
1207 bbox = geom.BoxI(
1208 geom.PointI(int(center[0] - 0.5 * width), int(center[1] - 0.5 * height)),
1209 geom.PointI(int(center[0] + 0.5 * width), int(center[1] + 0.5 * height)),
1210 )
1211 spans = afwGeom.SpanSet(bbox)
1212 elif rec["type"] == "circle":
1213 radius = int(rec["radius"].asArcseconds() / plateScale) # convert to pixels
1214 spans = afwGeom.SpanSet.fromShape(radius, offset=center)
1215 else:
1216 self.log.warning("Unexpected region type %s at %s", rec["type"], center)
1217 continue
1218 spans.clippedTo(mask.getBBox()).setMask(mask, self.brightObjectBitmask)
1220 def setInexactPsf(self, mask):
1221 """Set INEXACT_PSF mask plane.
1223 If any of the input images isn't represented in the coadd (due to
1224 clipped pixels or chip gaps), the `CoaddPsf` will be inexact. Flag
1225 these pixels.
1227 Parameters
1228 ----------
1229 mask : `lsst.afw.image.Mask`
1230 Coadded exposure's mask, modified in-place.
1231 """
1232 mask.addMaskPlane("INEXACT_PSF")
1233 inexactPsf = mask.getPlaneBitMask("INEXACT_PSF")
1234 sensorEdge = mask.getPlaneBitMask("SENSOR_EDGE") # chip edges (so PSF is discontinuous)
1235 clipped = mask.getPlaneBitMask("CLIPPED") # pixels clipped from coadd
1236 rejected = mask.getPlaneBitMask("REJECTED") # pixels rejected from coadd due to masks
1237 array = mask.getArray()
1238 selected = array & (sensorEdge | clipped | rejected) > 0
1239 array[selected] |= inexactPsf
1241 def filterWarps(self, inputs, goodVisits):
1242 """Return list of only inputRefs with visitId in goodVisits ordered by
1243 goodVisit.
1245 Parameters
1246 ----------
1247 inputs : `list` of `~lsst.pipe.base.connections.DeferredDatasetRef`
1248 List of `lsst.pipe.base.connections.DeferredDatasetRef` with dataId
1249 containing visit.
1250 goodVisit : `dict`
1251 Dictionary with good visitIds as the keys. Value ignored.
1253 Returns
1254 -------
1255 filterInputs : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]
1256 Filtered and sorted list of inputRefs with visitId in goodVisits
1257 ordered by goodVisit.
1258 """
1259 inputWarpDict = {inputRef.ref.dataId["visit"]: inputRef for inputRef in inputs}
1260 filteredInputs = []
1261 for visit in goodVisits.keys():
1262 if visit in inputWarpDict:
1263 filteredInputs.append(inputWarpDict[visit])
1264 return filteredInputs
1267def countMaskFromFootprint(mask, footprint, bitmask, ignoreMask):
1268 """Function to count the number of pixels with a specific mask in a
1269 footprint.
1271 Find the intersection of mask & footprint. Count all pixels in the mask
1272 that are in the intersection that have bitmask set but do not have
1273 ignoreMask set. Return the count.
1275 Parameters
1276 ----------
1277 mask : `lsst.afw.image.Mask`
1278 Mask to define intersection region by.
1279 footprint : `lsst.afw.detection.Footprint`
1280 Footprint to define the intersection region by.
1281 bitmask : `Unknown`
1282 Specific mask that we wish to count the number of occurances of.
1283 ignoreMask : `Unknown`
1284 Pixels to not consider.
1286 Returns
1287 -------
1288 result : `int`
1289 Number of pixels in footprint with specified mask.
1290 """
1291 bbox = footprint.getBBox()
1292 bbox.clip(mask.getBBox(afwImage.PARENT))
1293 fp = afwImage.Mask(bbox)
1294 subMask = mask.Factory(mask, bbox, afwImage.PARENT)
1295 footprint.spans.setMask(fp, bitmask)
1296 return numpy.logical_and(
1297 (subMask.getArray() & fp.getArray()) > 0, (subMask.getArray() & ignoreMask) == 0
1298 ).sum()
1301class CompareWarpAssembleCoaddConnections(AssembleCoaddConnections):
1302 psfMatchedWarps = pipeBase.connectionTypes.Input(
1303 doc=(
1304 "PSF-Matched Warps are required by CompareWarp regardless of the coadd type requested. "
1305 "Only PSF-Matched Warps make sense for image subtraction. "
1306 "Therefore, they must be an additional declared input."
1307 ),
1308 name="{inputCoaddName}Coadd_psfMatchedWarp",
1309 storageClass="ExposureF",
1310 dimensions=("tract", "patch", "skymap", "visit"),
1311 deferLoad=True,
1312 multiple=True,
1313 )
1314 templateCoadd = pipeBase.connectionTypes.Output(
1315 doc=(
1316 "Model of the static sky, used to find temporal artifacts. Typically a PSF-Matched, "
1317 "sigma-clipped coadd. Written if and only if assembleStaticSkyModel.doWrite=True"
1318 ),
1319 name="{outputCoaddName}CoaddPsfMatched",
1320 storageClass="ExposureF",
1321 dimensions=("tract", "patch", "skymap", "band"),
1322 )
1323 artifactMasks = pipeBase.connectionTypes.Output(
1324 doc="Mask of artifacts detected in the coadd",
1325 name="compare_warp_artifact_mask",
1326 storageClass="Mask",
1327 dimensions=("tract", "patch", "skymap", "visit", "instrument"),
1328 multiple=True,
1329 )
1331 def __init__(self, *, config=None):
1332 super().__init__(config=config)
1333 if not config.assembleStaticSkyModel.doWrite:
1334 self.outputs.remove("templateCoadd")
1335 config.validate()
1338class CompareWarpAssembleCoaddConfig(
1339 AssembleCoaddConfig, pipelineConnections=CompareWarpAssembleCoaddConnections
1340):
1341 assembleStaticSkyModel = pexConfig.ConfigurableField(
1342 target=AssembleCoaddTask,
1343 doc="Task to assemble an artifact-free, PSF-matched Coadd to serve as "
1344 "a naive/first-iteration model of the static sky.",
1345 )
1346 detect = pexConfig.ConfigurableField(
1347 target=SourceDetectionTask,
1348 doc="Detect outlier sources on difference between each psfMatched warp and static sky model",
1349 )
1350 detectTemplate = pexConfig.ConfigurableField(
1351 target=SourceDetectionTask,
1352 doc="Detect sources on static sky model. Only used if doPreserveContainedBySource is True",
1353 )
1354 maskStreaks = pexConfig.ConfigurableField(
1355 target=MaskStreaksTask,
1356 doc="Detect streaks on difference between each psfMatched warp and static sky model. Only used if "
1357 "doFilterMorphological is True. Adds a mask plane to an exposure, with the mask plane name set by"
1358 "streakMaskName",
1359 )
1360 streakMaskName = pexConfig.Field(dtype=str, default="STREAK", doc="Name of mask bit used for streaks")
1361 maxNumEpochs = pexConfig.Field(
1362 doc="Charactistic maximum local number of epochs/visits in which an artifact candidate can appear "
1363 "and still be masked. The effective maxNumEpochs is a broken linear function of local "
1364 "number of epochs (N): min(maxFractionEpochsLow*N, maxNumEpochs + maxFractionEpochsHigh*N). "
1365 "For each footprint detected on the image difference between the psfMatched warp and static sky "
1366 "model, if a significant fraction of pixels (defined by spatialThreshold) are residuals in more "
1367 "than the computed effective maxNumEpochs, the artifact candidate is deemed persistant rather "
1368 "than transient and not masked.",
1369 dtype=int,
1370 default=2,
1371 )
1372 maxFractionEpochsLow = pexConfig.RangeField(
1373 doc="Fraction of local number of epochs (N) to use as effective maxNumEpochs for low N. "
1374 "Effective maxNumEpochs = "
1375 "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)",
1376 dtype=float,
1377 default=0.4,
1378 min=0.0,
1379 max=1.0,
1380 )
1381 maxFractionEpochsHigh = pexConfig.RangeField(
1382 doc="Fraction of local number of epochs (N) to use as effective maxNumEpochs for high N. "
1383 "Effective maxNumEpochs = "
1384 "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)",
1385 dtype=float,
1386 default=0.03,
1387 min=0.0,
1388 max=1.0,
1389 )
1390 spatialThreshold = pexConfig.RangeField(
1391 doc="Unitless fraction of pixels defining how much of the outlier region has to meet the "
1392 "temporal criteria. If 0, clip all. If 1, clip none.",
1393 dtype=float,
1394 default=0.5,
1395 min=0.0,
1396 max=1.0,
1397 inclusiveMin=True,
1398 inclusiveMax=True,
1399 )
1400 doScaleWarpVariance = pexConfig.Field(
1401 doc="Rescale Warp variance plane using empirical noise?",
1402 dtype=bool,
1403 default=True,
1404 )
1405 scaleWarpVariance = pexConfig.ConfigurableField(
1406 target=ScaleVarianceTask,
1407 doc="Rescale variance on warps",
1408 )
1409 doPreserveContainedBySource = pexConfig.Field(
1410 doc="Rescue artifacts from clipping that completely lie within a footprint detected"
1411 "on the PsfMatched Template Coadd. Replicates a behavior of SafeClip.",
1412 dtype=bool,
1413 default=True,
1414 )
1415 doPrefilterArtifacts = pexConfig.Field(
1416 doc="Ignore artifact candidates that are mostly covered by the bad pixel mask, "
1417 "because they will be excluded anyway. This prevents them from contributing "
1418 "to the outlier epoch count image and potentially being labeled as persistant."
1419 "'Mostly' is defined by the config 'prefilterArtifactsRatio'.",
1420 dtype=bool,
1421 default=True,
1422 )
1423 prefilterArtifactsMaskPlanes = pexConfig.ListField(
1424 doc="Prefilter artifact candidates that are mostly covered by these bad mask planes.",
1425 dtype=str,
1426 default=("NO_DATA", "BAD", "SUSPECT"),
1427 )
1428 prefilterArtifactsRatio = pexConfig.Field(
1429 doc="Prefilter artifact candidates with less than this fraction overlapping good pixels",
1430 dtype=float,
1431 default=0.05,
1432 )
1433 doFilterMorphological = pexConfig.Field(
1434 doc="Filter artifact candidates based on morphological criteria, i.g. those that appear to "
1435 "be streaks.",
1436 dtype=bool,
1437 default=False,
1438 )
1439 growStreakFp = pexConfig.Field(
1440 doc="Grow streak footprints by this number multiplied by the PSF width", dtype=float, default=5
1441 )
1443 def setDefaults(self):
1444 AssembleCoaddConfig.setDefaults(self)
1445 self.statistic = "MEAN"
1447 # Real EDGE removed by psfMatched NO_DATA border half the width of the
1448 # matching kernel. CompareWarp applies psfMatched EDGE pixels to
1449 # directWarps before assembling.
1450 if "EDGE" in self.badMaskPlanes:
1451 self.badMaskPlanes.remove("EDGE")
1452 self.removeMaskPlanes.append("EDGE")
1453 self.assembleStaticSkyModel.badMaskPlanes = [
1454 "NO_DATA",
1455 ]
1456 self.assembleStaticSkyModel.warpType = "psfMatched"
1457 self.assembleStaticSkyModel.connections.warpType = "psfMatched"
1458 self.assembleStaticSkyModel.statistic = "MEANCLIP"
1459 self.assembleStaticSkyModel.sigmaClip = 2.5
1460 self.assembleStaticSkyModel.clipIter = 3
1461 self.assembleStaticSkyModel.calcErrorFromInputVariance = True
1462 self.assembleStaticSkyModel.doWrite = False
1463 self.detect.doTempLocalBackground = False
1464 self.detect.reEstimateBackground = False
1465 self.detect.returnOriginalFootprints = False
1466 self.detect.thresholdPolarity = "both"
1467 self.detect.thresholdValue = 7.5
1468 self.detect.minPixels = 4
1469 self.detect.isotropicGrow = True
1470 self.detect.thresholdType = "pixel_stdev"
1471 self.detect.nSigmaToGrow = 0.4
1472 # The default nSigmaToGrow for SourceDetectionTask is already 2.4,
1473 # Explicitly restating because ratio with detect.nSigmaToGrow matters
1474 self.detectTemplate.nSigmaToGrow = 2.4
1475 # Higher thresholds make smaller and fewer protected zones around
1476 # bright stars. Sources with snr < 300 tend to subtract OK empirically
1477 self.detectTemplate.thresholdValue = 300
1478 self.detectTemplate.doTempLocalBackground = True
1479 self.detectTemplate.reEstimateBackground = True
1480 self.detectTemplate.background.useApprox = False
1481 self.detectTemplate.returnOriginalFootprints = False
1483 def validate(self):
1484 super().validate()
1485 if self.assembleStaticSkyModel.doNImage:
1486 raise ValueError(
1487 "No dataset type exists for a PSF-Matched Template N Image."
1488 "Please set assembleStaticSkyModel.doNImage=False"
1489 )
1491 if self.assembleStaticSkyModel.doWrite and (self.warpType == self.assembleStaticSkyModel.warpType):
1492 raise ValueError(
1493 "warpType (%s) == assembleStaticSkyModel.warpType (%s) and will compete for "
1494 "the same dataset name. Please set assembleStaticSkyModel.doWrite to False "
1495 "or warpType to 'direct'. assembleStaticSkyModel.warpType should ways be "
1496 "'PsfMatched'" % (self.warpType, self.assembleStaticSkyModel.warpType)
1497 )
1500class CompareWarpAssembleCoaddTask(AssembleCoaddTask):
1501 """Assemble a compareWarp coadded image from a set of warps
1502 by masking artifacts detected by comparing PSF-matched warps.
1504 In ``AssembleCoaddTask``, we compute the coadd as an clipped mean (i.e.,
1505 we clip outliers). The problem with doing this is that when computing the
1506 coadd PSF at a given location, individual visit PSFs from visits with
1507 outlier pixels contribute to the coadd PSF and cannot be treated correctly.
1508 In this task, we correct for this behavior by creating a new badMaskPlane
1509 'CLIPPED' which marks pixels in the individual warps suspected to contain
1510 an artifact. We populate this plane on the input warps by comparing
1511 PSF-matched warps with a PSF-matched median coadd which serves as a
1512 model of the static sky. Any group of pixels that deviates from the
1513 PSF-matched template coadd by more than config.detect.threshold sigma,
1514 is an artifact candidate. The candidates are then filtered to remove
1515 variable sources and sources that are difficult to subtract such as
1516 bright stars. This filter is configured using the config parameters
1517 ``temporalThreshold`` and ``spatialThreshold``. The temporalThreshold is
1518 the maximum fraction of epochs that the deviation can appear in and still
1519 be considered an artifact. The spatialThreshold is the maximum fraction of
1520 pixels in the footprint of the deviation that appear in other epochs
1521 (where other epochs is defined by the temporalThreshold). If the deviant
1522 region meets this criteria of having a significant percentage of pixels
1523 that deviate in only a few epochs, these pixels have the 'CLIPPED' bit
1524 set in the mask. These regions will not contribute to the final coadd.
1525 Furthermore, any routine to determine the coadd PSF can now be cognizant
1526 of clipped regions. Note that the algorithm implemented by this task is
1527 preliminary and works correctly for HSC data. Parameter modifications and
1528 or considerable redesigning of the algorithm is likley required for other
1529 surveys.
1531 ``CompareWarpAssembleCoaddTask`` sub-classes
1532 ``AssembleCoaddTask`` and instantiates ``AssembleCoaddTask``
1533 as a subtask to generate the TemplateCoadd (the model of the static sky).
1535 Notes
1536 -----
1537 Debugging:
1538 This task supports the following debug variables:
1539 - ``saveCountIm``
1540 If True then save the Epoch Count Image as a fits file in the `figPath`
1541 - ``figPath``
1542 Path to save the debug fits images and figures
1543 """
1545 ConfigClass = CompareWarpAssembleCoaddConfig
1546 _DefaultName = "compareWarpAssembleCoadd"
1548 # See the parent class for docstring.
1549 _doUsePsfMatchedPolygons: bool = True
1551 def __init__(self, *args, **kwargs):
1552 AssembleCoaddTask.__init__(self, *args, **kwargs)
1553 self.makeSubtask("assembleStaticSkyModel")
1554 detectionSchema = afwTable.SourceTable.makeMinimalSchema()
1555 self.makeSubtask("detect", schema=detectionSchema)
1556 if self.config.doPreserveContainedBySource:
1557 self.makeSubtask("detectTemplate", schema=afwTable.SourceTable.makeMinimalSchema())
1558 if self.config.doScaleWarpVariance:
1559 self.makeSubtask("scaleWarpVariance")
1560 if self.config.doFilterMorphological:
1561 self.makeSubtask("maskStreaks")
1563 @utils.inheritDoc(AssembleCoaddTask)
1564 def _makeSupplementaryData(self, butlerQC, inputRefs, outputRefs):
1565 """Generate a templateCoadd to use as a naive model of static sky to
1566 subtract from PSF-Matched warps.
1568 Returns
1569 -------
1570 result : `~lsst.pipe.base.Struct`
1571 Results as a struct with attributes:
1573 ``templateCoadd``
1574 Coadded exposure (`lsst.afw.image.Exposure`).
1575 ``nImage``
1576 Keeps track of exposure count for each pixel
1577 (`lsst.afw.image.ImageU`).
1579 Raises
1580 ------
1581 RuntimeError
1582 Raised if ``templateCoadd`` is `None`.
1583 """
1584 # Ensure that psfMatchedWarps are used as input warps for template
1585 # generation.
1586 staticSkyModelInputRefs = copy.deepcopy(inputRefs)
1587 staticSkyModelInputRefs.inputWarps = inputRefs.psfMatchedWarps
1589 # self.assembleStaticSkyModel is not expecting psfMatchedWarps as
1590 # input. But in its runQuantum, ButlerQC would try to be helpful and
1591 # read all of them in to memory.
1592 del staticSkyModelInputRefs.psfMatchedWarps
1594 # Because subtasks don't have connections we have to make one.
1595 # The main task's `templateCoadd` is the subtask's `coaddExposure`
1596 staticSkyModelOutputRefs = copy.deepcopy(outputRefs)
1597 if self.config.assembleStaticSkyModel.doWrite:
1598 staticSkyModelOutputRefs.coaddExposure = staticSkyModelOutputRefs.templateCoadd
1599 # Remove template coadd from both subtask's and main tasks outputs,
1600 # because it is handled by the subtask as `coaddExposure`
1601 del outputRefs.templateCoadd
1602 del staticSkyModelOutputRefs.templateCoadd
1604 # A PSF-Matched nImage does not exist as a dataset type
1605 if "nImage" in staticSkyModelOutputRefs.keys():
1606 del staticSkyModelOutputRefs.nImage
1608 templateCoadd = self.assembleStaticSkyModel.runQuantum(
1609 butlerQC, staticSkyModelInputRefs, staticSkyModelOutputRefs
1610 )
1611 if templateCoadd is None:
1612 raise RuntimeError(self._noTemplateMessage(self.assembleStaticSkyModel.warpType))
1614 return pipeBase.Struct(
1615 templateCoadd=templateCoadd.coaddExposure,
1616 nImage=templateCoadd.nImage,
1617 warpRefList=templateCoadd.warpRefList,
1618 imageScalerList=templateCoadd.imageScalerList,
1619 weightList=templateCoadd.weightList,
1620 psfMatchedWarpRefList=inputRefs.psfMatchedWarps,
1621 )
1623 def _noTemplateMessage(self, warpType):
1624 warpName = warpType[0].upper() + warpType[1:]
1625 message = """No %(warpName)s warps were found to build the template coadd which is
1626 required to run CompareWarpAssembleCoaddTask. To continue assembling this type of coadd,
1627 first either rerun makeCoaddTempExp with config.make%(warpName)s=True or
1628 coaddDriver with config.makeCoadTempExp.make%(warpName)s=True, before assembleCoadd.
1630 Alternatively, to use another algorithm with existing warps, retarget the CoaddDriverConfig to
1631 another algorithm like:
1633 from lsst.pipe.tasks.assembleCoadd import SafeClipAssembleCoaddTask
1634 config.assemble.retarget(SafeClipAssembleCoaddTask)
1635 """ % {
1636 "warpName": warpName
1637 }
1638 return message
1640 @utils.inheritDoc(AssembleCoaddTask)
1641 @timeMethod
1642 def run(
1643 self,
1644 skyInfo,
1645 *,
1646 warpRefList,
1647 imageScalerList,
1648 weightList,
1649 psfMatchedWarpRefList,
1650 supplementaryData,
1651 ):
1652 """Notes
1653 -----
1654 Assemble the coadd.
1656 Find artifacts and apply them to the warps' masks creating a list of
1657 alternative masks with a new "CLIPPED" plane and updated "NO_DATA"
1658 plane. Then pass these alternative masks to the base class's ``run``
1659 method.
1660 """
1661 # Check and match the order of the supplementaryData
1662 # (PSF-matched) inputs to the order of the direct inputs,
1663 # so that the artifact mask is applied to the right warp
1665 # For any missing psfMatched warp refs replaced with None,
1666 # findArtifacts() will produce a NO_DATA mask covering the entire bbox.
1667 # As long as NO_DATA is in self.config.badMaskPlanes, these direct
1668 # warps with missing psfMatched warps will not contribute to the coadd.
1669 dataIds = [ref.dataId for ref in warpRefList]
1670 psfMatchedDataIds = [ref.dataId for ref in supplementaryData.warpRefList]
1672 if dataIds != psfMatchedDataIds:
1673 self.log.info("Reordering and or/padding PSF-matched visit input list")
1674 supplementaryData.warpRefList = reorderAndPadList(
1675 supplementaryData.warpRefList, psfMatchedDataIds, dataIds
1676 )
1677 # Remove in DM-49083
1678 supplementaryData.imageScalerList = reorderAndPadList(
1679 supplementaryData.imageScalerList, psfMatchedDataIds, dataIds
1680 )
1682 # Use PSF-Matched Warps (and corresponding scalers) and coadd to find
1683 # artifacts.
1684 spanSetMaskList = self.findArtifacts(
1685 supplementaryData.templateCoadd, supplementaryData.warpRefList, supplementaryData.imageScalerList
1686 )
1687 # The mask planes are defined globally, so we can load the dict from
1688 # the templateCoadd or the warps.
1689 badMaskPlanes = [
1690 mp
1691 for mp in self.config.badMaskPlanes
1692 if mp in supplementaryData.templateCoadd.mask.getMaskPlaneDict().keys()
1693 ]
1694 badMaskPlanes.append("CLIPPED")
1695 badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes)
1697 result = AssembleCoaddTask.run(
1698 self,
1699 skyInfo,
1700 warpRefList=warpRefList,
1701 imageScalerList=imageScalerList,
1702 weightList=weightList,
1703 altMaskList=spanSetMaskList,
1704 mask=badPixelMask,
1705 psfMatchedWarpRefList=psfMatchedWarpRefList,
1706 )
1708 # Propagate PSF-matched EDGE pixels to coadd SENSOR_EDGE and
1709 # INEXACT_PSF. Psf-Matching moves the real edge inwards.
1710 self.applyAltEdgeMask(result.coaddExposure.maskedImage.mask, spanSetMaskList)
1711 return result
1713 def applyAltEdgeMask(self, mask, altMaskList):
1714 """Propagate alt EDGE mask to SENSOR_EDGE AND INEXACT_PSF planes.
1716 Parameters
1717 ----------
1718 mask : `lsst.afw.image.Mask`
1719 Original mask.
1720 altMaskList : `list` of `dict`
1721 List of Dicts containing ``spanSet`` lists.
1722 Each element contains the new mask plane name (e.g. "CLIPPED
1723 and/or "NO_DATA") as the key, and list of ``SpanSets`` to apply to
1724 the mask.
1725 """
1726 maskValue = mask.getPlaneBitMask(["SENSOR_EDGE", "INEXACT_PSF"])
1727 for visitMask in altMaskList:
1728 if "EDGE" in visitMask:
1729 for spanSet in visitMask["EDGE"]:
1730 spanSet.clippedTo(mask.getBBox()).setMask(mask, maskValue)
1732 def findArtifacts(self, templateCoadd, warpRefList, imageScalerList):
1733 """Find artifacts.
1735 Loop through warps twice. The first loop builds a map with the count
1736 of how many epochs each pixel deviates from the templateCoadd by more
1737 than ``config.chiThreshold`` sigma. The second loop takes each
1738 difference image and filters the artifacts detected in each using
1739 count map to filter out variable sources and sources that are
1740 difficult to subtract cleanly.
1742 Parameters
1743 ----------
1744 templateCoadd : `lsst.afw.image.Exposure`
1745 Exposure to serve as model of static sky.
1746 warpRefList : `list`
1747 List of dataset handles (data references) to warps.
1748 imageScalerList : `list`
1749 List of image scalers.
1750 Deprecated and will be removed after v29.
1752 Returns
1753 -------
1754 altMasks : `list` of `dict`
1755 List of dicts containing information about CLIPPED
1756 (i.e., artifacts), NO_DATA, and EDGE pixels.
1757 """
1758 self.log.debug("Generating Count Image, and mask lists.")
1759 coaddBBox = templateCoadd.getBBox()
1760 slateIm = afwImage.ImageU(coaddBBox)
1761 epochCountImage = afwImage.ImageU(coaddBBox)
1762 nImage = afwImage.ImageU(coaddBBox)
1763 spanSetArtifactList = []
1764 spanSetNoDataMaskList = []
1765 spanSetEdgeList = []
1766 spanSetBadMorphoList = []
1767 badPixelMask = self.getBadPixelMask()
1769 # mask of the warp diffs should = that of only the warp
1770 templateCoadd.mask.clearAllMaskPlanes()
1772 if self.config.doPreserveContainedBySource:
1773 sacrificeToDetection = templateCoadd.clone()
1774 templateFootprints = self.detectTemplate.detectFootprints(sacrificeToDetection)
1775 del sacrificeToDetection
1776 else:
1777 templateFootprints = None
1779 for warpRef, imageScaler in zip(warpRefList, imageScalerList):
1780 warpDiffExp = self._readAndComputeWarpDiff(warpRef, imageScaler, templateCoadd)
1781 if warpDiffExp is not None:
1782 # This nImage only approximates the final nImage because it
1783 # uses the PSF-matched mask.
1784 nImage.array += (
1785 numpy.isfinite(warpDiffExp.image.array) * ((warpDiffExp.mask.array & badPixelMask) == 0)
1786 ).astype(numpy.uint16)
1787 fpSet = self.detect.detectFootprints(warpDiffExp, doSmooth=False, clearMask=True)
1788 fpSet.positive.merge(fpSet.negative)
1789 footprints = fpSet.positive
1790 slateIm.set(0)
1791 spanSetList = [footprint.spans for footprint in footprints.getFootprints()]
1793 # Remove artifacts due to defects before they contribute to
1794 # the epochCountImage.
1795 if self.config.doPrefilterArtifacts:
1796 spanSetList = self.prefilterArtifacts(spanSetList, warpDiffExp)
1798 # Clear mask before adding prefiltered spanSets
1799 self.detect.clearMask(warpDiffExp.mask)
1800 for spans in spanSetList:
1801 spans.setImage(slateIm, 1, doClip=True)
1802 spans.setMask(warpDiffExp.mask, warpDiffExp.mask.getPlaneBitMask("DETECTED"))
1803 epochCountImage += slateIm
1805 if self.config.doFilterMorphological:
1806 maskName = self.config.streakMaskName
1807 # clear single frame streak mask if it exists already
1808 if maskName in warpDiffExp.mask.getMaskPlaneDict():
1809 warpDiffExp.mask.clearMaskPlane(warpDiffExp.mask.getMaskPlane(maskName))
1810 else:
1811 self.log.debug(f"Did not (need to) clear {maskName} mask because it didn't exist")
1813 _ = self.maskStreaks.run(warpDiffExp)
1814 streakMask = warpDiffExp.mask
1815 initSpanSetStreak = afwGeom.SpanSet.fromMask(
1816 streakMask, streakMask.getPlaneBitMask(maskName)
1817 ).split()
1818 # Pad the streaks to account for low-surface brightness
1819 # wings.
1820 psf = warpDiffExp.getPsf()
1821 spanSetStreak = []
1822 for sset in initSpanSetStreak:
1823 if self.config.doPreserveContainedBySource and templateFootprints is not None:
1824 doKeep = True
1825 for footprint in templateFootprints.positive.getFootprints():
1826 if footprint.spans.contains(sset):
1827 doKeep = False
1828 break
1829 if not doKeep:
1830 continue
1831 psfShape = psf.computeShape(sset.computeCentroid())
1832 dilation = self.config.growStreakFp * psfShape.getDeterminantRadius()
1833 sset_dilated = sset.dilated(int(dilation))
1834 spanSetStreak.append(sset_dilated)
1836 # PSF-Matched warps have less available area (~the matching
1837 # kernel) because the calexps undergo a second convolution.
1838 # Pixels with data in the direct warp but not in the
1839 # PSF-matched warp will not have their artifacts detected.
1840 # NaNs from the PSF-matched warp therefore must be masked in
1841 # the direct warp.
1842 nans = numpy.where(numpy.isnan(warpDiffExp.maskedImage.image.array), 1, 0)
1843 nansMask = afwImage.makeMaskFromArray(nans.astype(afwImage.MaskPixel))
1844 nansMask.setXY0(warpDiffExp.getXY0())
1845 edgeMask = warpDiffExp.mask
1846 spanSetEdgeMask = afwGeom.SpanSet.fromMask(edgeMask, edgeMask.getPlaneBitMask("EDGE")).split()
1847 else:
1848 # If the directWarp has <1% coverage, the psfMatchedWarp can
1849 # have 0% and not exist. In this case, mask the whole epoch.
1850 nansMask = afwImage.MaskX(coaddBBox, 1)
1851 spanSetList = []
1852 spanSetEdgeMask = []
1853 spanSetStreak = []
1855 spanSetNoDataMask = afwGeom.SpanSet.fromMask(nansMask).split()
1857 spanSetNoDataMaskList.append(spanSetNoDataMask)
1858 spanSetArtifactList.append(spanSetList)
1859 spanSetEdgeList.append(spanSetEdgeMask)
1860 if self.config.doFilterMorphological:
1861 spanSetBadMorphoList.append(spanSetStreak)
1863 if lsstDebug.Info(__name__).saveCountIm:
1864 path = self._dataRef2DebugPath("epochCountIm", warpRefList[0], coaddLevel=True)
1865 epochCountImage.writeFits(path)
1867 for i, spanSetList in enumerate(spanSetArtifactList):
1868 if spanSetList:
1869 filteredSpanSetList = self.filterArtifacts(
1870 spanSetList, epochCountImage, nImage, templateFootprints
1871 )
1872 spanSetArtifactList[i] = filteredSpanSetList
1873 if self.config.doFilterMorphological:
1874 spanSetArtifactList[i] += spanSetBadMorphoList[i]
1876 altMasks = []
1877 for artifacts, noData, edge in zip(spanSetArtifactList, spanSetNoDataMaskList, spanSetEdgeList):
1878 altMasks.append({"CLIPPED": artifacts, "NO_DATA": noData, "EDGE": edge})
1879 return altMasks
1881 def prefilterArtifacts(self, spanSetList, exp):
1882 """Remove artifact candidates covered by bad mask plane.
1884 Any future editing of the candidate list that does not depend on
1885 temporal information should go in this method.
1887 Parameters
1888 ----------
1889 spanSetList : `list` [`lsst.afw.geom.SpanSet`]
1890 List of SpanSets representing artifact candidates.
1891 exp : `lsst.afw.image.Exposure`
1892 Exposure containing mask planes used to prefilter.
1894 Returns
1895 -------
1896 returnSpanSetList : `list` [`lsst.afw.geom.SpanSet`]
1897 List of SpanSets with artifacts.
1898 """
1899 existingMaskPlanes = [
1900 m for m in self.config.prefilterArtifactsMaskPlanes if m in exp.mask.getMaskPlaneDict()
1901 ]
1902 badPixelMask = exp.mask.getPlaneBitMask(existingMaskPlanes)
1903 goodArr = (exp.mask.array & badPixelMask) == 0
1904 returnSpanSetList = []
1905 bbox = exp.getBBox()
1906 x0, y0 = exp.getXY0()
1907 for i, span in enumerate(spanSetList):
1908 y, x = span.clippedTo(bbox).indices()
1909 yIndexLocal = numpy.array(y) - y0
1910 xIndexLocal = numpy.array(x) - x0
1911 goodRatio = numpy.count_nonzero(goodArr[yIndexLocal, xIndexLocal]) / span.getArea()
1912 if goodRatio > self.config.prefilterArtifactsRatio:
1913 returnSpanSetList.append(span)
1914 return returnSpanSetList
1916 def filterArtifacts(self, spanSetList, epochCountImage, nImage, footprintsToExclude=None):
1917 """Filter artifact candidates.
1919 Parameters
1920 ----------
1921 spanSetList : `list` [`lsst.afw.geom.SpanSet`]
1922 List of SpanSets representing artifact candidates.
1923 epochCountImage : `lsst.afw.image.Image`
1924 Image of accumulated number of warpDiff detections.
1925 nImage : `lsst.afw.image.ImageU`
1926 Image of the accumulated number of total epochs contributing.
1928 Returns
1929 -------
1930 maskSpanSetList : `list` [`lsst.afw.geom.SpanSet`]
1931 List of SpanSets with artifacts.
1932 """
1933 maskSpanSetList = []
1934 x0, y0 = epochCountImage.getXY0()
1935 for i, span in enumerate(spanSetList):
1936 y, x = span.indices()
1937 yIdxLocal = [y1 - y0 for y1 in y]
1938 xIdxLocal = [x1 - x0 for x1 in x]
1939 outlierN = epochCountImage.array[yIdxLocal, xIdxLocal]
1940 totalN = nImage.array[yIdxLocal, xIdxLocal]
1942 # effectiveMaxNumEpochs is broken line (fraction of N) with
1943 # characteristic config.maxNumEpochs.
1944 effMaxNumEpochsHighN = self.config.maxNumEpochs + self.config.maxFractionEpochsHigh * numpy.mean(
1945 totalN
1946 )
1947 effMaxNumEpochsLowN = self.config.maxFractionEpochsLow * numpy.mean(totalN)
1948 effectiveMaxNumEpochs = int(min(effMaxNumEpochsLowN, effMaxNumEpochsHighN))
1949 nPixelsBelowThreshold = numpy.count_nonzero((outlierN > 0) & (outlierN <= effectiveMaxNumEpochs))
1950 percentBelowThreshold = nPixelsBelowThreshold / len(outlierN)
1951 if percentBelowThreshold > self.config.spatialThreshold:
1952 maskSpanSetList.append(span)
1954 if self.config.doPreserveContainedBySource and footprintsToExclude is not None:
1955 # If a candidate is contained by a footprint on the template coadd,
1956 # do not clip.
1957 filteredMaskSpanSetList = []
1958 for span in maskSpanSetList:
1959 doKeep = True
1960 for footprint in footprintsToExclude.positive.getFootprints():
1961 if footprint.spans.contains(span):
1962 doKeep = False
1963 break
1964 if doKeep:
1965 filteredMaskSpanSetList.append(span)
1966 maskSpanSetList = filteredMaskSpanSetList
1968 return maskSpanSetList
1970 def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd):
1971 """Fetch a warp from the butler and return a warpDiff.
1973 Parameters
1974 ----------
1975 warpRef : `lsst.daf.butler.DeferredDatasetHandle`
1976 Dataset handle for the warp.
1977 imageScaler : `lsst.pipe.tasks.scaleZeroPoint.ImageScaler`
1978 An image scaler object.
1979 Deprecated and will be removed after v29 in DM-49083.
1980 templateCoadd : `lsst.afw.image.Exposure`
1981 Exposure to be substracted from the scaled warp.
1983 Returns
1984 -------
1985 warp : `lsst.afw.image.Exposure`
1986 Exposure of the image difference between the warp and template.
1987 """
1988 # If the PSF-Matched warp did not exist for this direct warp
1989 # None is holding its place to maintain order in Gen 3
1990 if warpRef is None:
1991 return None
1993 warp = warpRef.get(parameters={"bbox": templateCoadd.getBBox()})
1994 # direct image scaler OK for PSF-matched Warp.
1995 # Remove in DM-49083.
1996 if imageScaler is not None:
1997 imageScaler.scaleMaskedImage(warp.getMaskedImage())
1998 mi = warp.getMaskedImage()
1999 if self.config.doScaleWarpVariance:
2000 try:
2001 self.scaleWarpVariance.run(mi)
2002 except Exception as exc:
2003 self.log.warning("Unable to rescale variance of warp (%s); leaving it as-is", exc)
2004 mi -= templateCoadd.getMaskedImage()
2005 return warp