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