41 from .coaddBase
import CoaddBaseTask, SelectDataIdContainer, makeSkyInfo, makeCoaddSuffix, reorderAndPadList
42 from .interpImage
import InterpImageTask
43 from .scaleZeroPoint
import ScaleZeroPointTask
44 from .coaddHelpers
import groupPatchExposures, getGroupDataRef
45 from .scaleVariance
import ScaleVarianceTask
46 from .maskStreaks
import MaskStreaksTask
47 from .healSparseMapping
import HealSparseInputMapTask
49 from lsst.daf.butler
import DeferredDatasetHandle
50 from lsst.utils.timer
import timeMethod
52 __all__ = [
"AssembleCoaddTask",
"AssembleCoaddConnections",
"AssembleCoaddConfig",
53 "SafeClipAssembleCoaddTask",
"SafeClipAssembleCoaddConfig",
54 "CompareWarpAssembleCoaddTask",
"CompareWarpAssembleCoaddConfig"]
56 log = logging.getLogger(__name__)
60 dimensions=(
"tract",
"patch",
"band",
"skymap"),
61 defaultTemplates={
"inputCoaddName":
"deep",
62 "outputCoaddName":
"deep",
64 "warpTypeSuffix":
""}):
66 inputWarps = pipeBase.connectionTypes.Input(
67 doc=(
"Input list of warps to be assemebled i.e. stacked."
68 "WarpType (e.g. direct, psfMatched) is controlled by the warpType config parameter"),
69 name=
"{inputCoaddName}Coadd_{warpType}Warp",
70 storageClass=
"ExposureF",
71 dimensions=(
"tract",
"patch",
"skymap",
"visit",
"instrument"),
75 skyMap = pipeBase.connectionTypes.Input(
76 doc=
"Input definition of geometry/bbox and projection/wcs for coadded exposures",
77 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
78 storageClass=
"SkyMap",
79 dimensions=(
"skymap", ),
81 selectedVisits = pipeBase.connectionTypes.Input(
82 doc=
"Selected visits to be coadded.",
83 name=
"{outputCoaddName}Visits",
84 storageClass=
"StructuredDataDict",
85 dimensions=(
"instrument",
"tract",
"patch",
"skymap",
"band")
87 brightObjectMask = pipeBase.connectionTypes.PrerequisiteInput(
88 doc=(
"Input Bright Object Mask mask produced with external catalogs to be applied to the mask plane"
90 name=
"brightObjectMask",
91 storageClass=
"ObjectMaskCatalog",
92 dimensions=(
"tract",
"patch",
"skymap",
"band"),
94 coaddExposure = pipeBase.connectionTypes.Output(
95 doc=
"Output coadded exposure, produced by stacking input warps",
96 name=
"{outputCoaddName}Coadd{warpTypeSuffix}",
97 storageClass=
"ExposureF",
98 dimensions=(
"tract",
"patch",
"skymap",
"band"),
100 nImage = pipeBase.connectionTypes.Output(
101 doc=
"Output image of number of input images per pixel",
102 name=
"{outputCoaddName}Coadd_nImage",
103 storageClass=
"ImageU",
104 dimensions=(
"tract",
"patch",
"skymap",
"band"),
106 inputMap = pipeBase.connectionTypes.Output(
107 doc=
"Output healsparse map of input images",
108 name=
"{outputCoaddName}Coadd_inputMap",
109 storageClass=
"HealSparseMap",
110 dimensions=(
"tract",
"patch",
"skymap",
"band"),
113 def __init__(self, *, config=None):
114 super().__init__(config=config)
119 templateValues = {name: getattr(config.connections, name)
for name
in self.defaultTemplates}
120 templateValues[
'warpType'] = config.warpType
122 self._nameOverrides = {name: getattr(config.connections, name).
format(**templateValues)
123 for name
in self.allConnections}
124 self._typeNameToVarName = {v: k
for k, v
in self._nameOverrides.
items()}
127 if not config.doMaskBrightObjects:
128 self.prerequisiteInputs.remove(
"brightObjectMask")
130 if not config.doSelectVisits:
131 self.inputs.remove(
"selectedVisits")
133 if not config.doNImage:
134 self.outputs.remove(
"nImage")
136 if not self.config.doInputMap:
137 self.outputs.remove(
"inputMap")
140 class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig,
141 pipelineConnections=AssembleCoaddConnections):
142 """Configuration parameters for the `AssembleCoaddTask`.
146 The `doMaskBrightObjects` and `brightObjectMaskName` configuration options
147 only set the bitplane config.brightObjectMaskName. To make this useful you
148 *must* also configure the flags.pixel algorithm, for example by adding
152 config.measurement.plugins["base_PixelFlags"].masksFpCenter.append("BRIGHT_OBJECT")
153 config.measurement.plugins["base_PixelFlags"].masksFpAnywhere.append("BRIGHT_OBJECT")
155 to your measureCoaddSources.py and forcedPhotCoadd.py config overrides.
157 warpType = pexConfig.Field(
158 doc=
"Warp name: one of 'direct' or 'psfMatched'",
162 subregionSize = pexConfig.ListField(
164 doc=
"Width, height of stack subregion size; "
165 "make small enough that a full stack of images will fit into memory at once.",
167 default=(2000, 2000),
169 statistic = pexConfig.Field(
171 doc=
"Main stacking statistic for aggregating over the epochs.",
174 doOnlineForMean = pexConfig.Field(
176 doc=
"Perform online coaddition when statistic=\"MEAN\" to save memory?",
179 doSigmaClip = pexConfig.Field(
181 doc=
"Perform sigma clipped outlier rejection with MEANCLIP statistic? (DEPRECATED)",
184 sigmaClip = pexConfig.Field(
186 doc=
"Sigma for outlier rejection; ignored if non-clipping statistic selected.",
189 clipIter = pexConfig.Field(
191 doc=
"Number of iterations of outlier rejection; ignored if non-clipping statistic selected.",
194 calcErrorFromInputVariance = pexConfig.Field(
196 doc=
"Calculate coadd variance from input variance by stacking statistic."
197 "Passed to StatisticsControl.setCalcErrorFromInputVariance()",
200 scaleZeroPoint = pexConfig.ConfigurableField(
201 target=ScaleZeroPointTask,
202 doc=
"Task to adjust the photometric zero point of the coadd temp exposures",
204 doInterp = pexConfig.Field(
205 doc=
"Interpolate over NaN pixels? Also extrapolate, if necessary, but the results are ugly.",
209 interpImage = pexConfig.ConfigurableField(
210 target=InterpImageTask,
211 doc=
"Task to interpolate (and extrapolate) over NaN pixels",
213 doWrite = pexConfig.Field(
214 doc=
"Persist coadd?",
218 doNImage = pexConfig.Field(
219 doc=
"Create image of number of contributing exposures for each pixel",
223 doUsePsfMatchedPolygons = pexConfig.Field(
224 doc=
"Use ValidPolygons from shrunk Psf-Matched Calexps? Should be set to True by CompareWarp only.",
228 maskPropagationThresholds = pexConfig.DictField(
231 doc=(
"Threshold (in fractional weight) of rejection at which we propagate a mask plane to "
232 "the coadd; that is, we set the mask bit on the coadd if the fraction the rejected frames "
233 "would have contributed exceeds this value."),
234 default={
"SAT": 0.1},
236 removeMaskPlanes = pexConfig.ListField(dtype=str, default=[
"NOT_DEBLENDED"],
237 doc=
"Mask planes to remove before coadding")
238 doMaskBrightObjects = pexConfig.Field(dtype=bool, default=
False,
239 doc=
"Set mask and flag bits for bright objects?")
240 brightObjectMaskName = pexConfig.Field(dtype=str, default=
"BRIGHT_OBJECT",
241 doc=
"Name of mask bit used for bright objects")
242 coaddPsf = pexConfig.ConfigField(
243 doc=
"Configuration for CoaddPsf",
244 dtype=measAlg.CoaddPsfConfig,
246 doAttachTransmissionCurve = pexConfig.Field(
247 dtype=bool, default=
False, optional=
False,
248 doc=(
"Attach a piecewise TransmissionCurve for the coadd? "
249 "(requires all input Exposures to have TransmissionCurves).")
251 hasFakes = pexConfig.Field(
254 doc=
"Should be set to True if fake sources have been inserted into the input data."
256 doSelectVisits = pexConfig.Field(
257 doc=
"Coadd only visits selected by a SelectVisitsTask",
261 doInputMap = pexConfig.Field(
262 doc=
"Create a bitwise map of coadd inputs",
266 inputMapper = pexConfig.ConfigurableField(
267 doc=
"Input map creation subtask.",
268 target=HealSparseInputMapTask,
273 self.badMaskPlanes = [
"NO_DATA",
"BAD",
"SAT",
"EDGE"]
280 log.warning(
"Config doPsfMatch deprecated. Setting warpType='psfMatched'")
281 self.warpType =
'psfMatched'
282 if self.doSigmaClip
and self.statistic !=
"MEANCLIP":
283 log.warning(
'doSigmaClip deprecated. To replicate behavior, setting statistic to "MEANCLIP"')
284 self.statistic =
"MEANCLIP"
285 if self.doInterp
and self.statistic
not in [
'MEAN',
'MEDIAN',
'MEANCLIP',
'VARIANCE',
'VARIANCECLIP']:
286 raise ValueError(
"Must set doInterp=False for statistic=%s, which does not "
287 "compute and set a non-zero coadd variance estimate." % (self.statistic))
289 unstackableStats = [
'NOTHING',
'ERROR',
'ORMASK']
290 if not hasattr(afwMath.Property, self.statistic)
or self.statistic
in unstackableStats:
291 stackableStats = [str(k)
for k
in afwMath.Property.__members__.keys()
292 if str(k)
not in unstackableStats]
293 raise ValueError(
"statistic %s is not allowed. Please choose one of %s."
294 % (self.statistic, stackableStats))
297 class AssembleCoaddTask(
CoaddBaseTask, pipeBase.PipelineTask):
298 """Assemble a coadded image from a set of warps (coadded temporary exposures).
300 We want to assemble a coadded image from a set of Warps (also called
301 coadded temporary exposures or ``coaddTempExps``).
302 Each input Warp covers a patch on the sky and corresponds to a single
303 run/visit/exposure of the covered patch. We provide the task with a list
304 of Warps (``selectDataList``) from which it selects Warps that cover the
305 specified patch (pointed at by ``dataRef``).
306 Each Warp that goes into a coadd will typically have an independent
307 photometric zero-point. Therefore, we must scale each Warp to set it to
308 a common photometric zeropoint. WarpType may be one of 'direct' or
309 'psfMatched', and the boolean configs `config.makeDirect` and
310 `config.makePsfMatched` set which of the warp types will be coadded.
311 The coadd is computed as a mean with optional outlier rejection.
312 Criteria for outlier rejection are set in `AssembleCoaddConfig`.
313 Finally, Warps can have bad 'NaN' pixels which received no input from the
314 source calExps. We interpolate over these bad (NaN) pixels.
316 `AssembleCoaddTask` uses several sub-tasks. These are
318 - `ScaleZeroPointTask`
319 - create and use an ``imageScaler`` object to scale the photometric zeropoint for each Warp
321 - interpolate across bad pixels (NaN) in the final coadd
323 You can retarget these subtasks if you wish.
327 The `lsst.pipe.base.cmdLineTask.CmdLineTask` interface supports a
328 flag ``-d`` to import ``debug.py`` from your ``PYTHONPATH``; see
329 `baseDebug` for more about ``debug.py`` files. `AssembleCoaddTask` has
330 no debug variables of its own. Some of the subtasks may support debug
331 variables. See the documentation for the subtasks for further information.
335 `AssembleCoaddTask` assembles a set of warped images into a coadded image.
336 The `AssembleCoaddTask` can be invoked by running ``assembleCoadd.py``
337 with the flag '--legacyCoadd'. Usage of assembleCoadd.py expects two
338 inputs: a data reference to the tract patch and filter to be coadded, and
339 a list of Warps to attempt to coadd. These are specified using ``--id`` and
340 ``--selectId``, respectively:
344 --id = [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]
345 --selectId [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]
347 Only the Warps that cover the specified tract and patch will be coadded.
348 A list of the available optional arguments can be obtained by calling
349 ``assembleCoadd.py`` with the ``--help`` command line argument:
353 assembleCoadd.py --help
355 To demonstrate usage of the `AssembleCoaddTask` in the larger context of
356 multi-band processing, we will generate the HSC-I & -R band coadds from
357 HSC engineering test data provided in the ``ci_hsc`` package. To begin,
358 assuming that the lsst stack has been already set up, we must set up the
359 obs_subaru and ``ci_hsc`` packages. This defines the environment variable
360 ``$CI_HSC_DIR`` and points at the location of the package. The raw HSC
361 data live in the ``$CI_HSC_DIR/raw directory``. To begin assembling the
362 coadds, we must first
365 - process the individual ccds in $CI_HSC_RAW to produce calibrated exposures
367 - create a skymap that covers the area of the sky present in the raw exposures
369 - warp the individual calibrated exposures to the tangent plane of the coadd
371 We can perform all of these steps by running
375 $CI_HSC_DIR scons warp-903986 warp-904014 warp-903990 warp-904010 warp-903988
377 This will produce warped exposures for each visit. To coadd the warped
378 data, we call assembleCoadd.py as follows:
382 assembleCoadd.py --legacyCoadd $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-I \
383 --selectId visit=903986 ccd=16 --selectId visit=903986 ccd=22 --selectId visit=903986 ccd=23 \
384 --selectId visit=903986 ccd=100 --selectId visit=904014 ccd=1 --selectId visit=904014 ccd=6 \
385 --selectId visit=904014 ccd=12 --selectId visit=903990 ccd=18 --selectId visit=903990 ccd=25 \
386 --selectId visit=904010 ccd=4 --selectId visit=904010 ccd=10 --selectId visit=904010 ccd=100 \
387 --selectId visit=903988 ccd=16 --selectId visit=903988 ccd=17 --selectId visit=903988 ccd=23 \
388 --selectId visit=903988 ccd=24
390 that will process the HSC-I band data. The results are written in
391 ``$CI_HSC_DIR/DATA/deepCoadd-results/HSC-I``.
393 You may also choose to run:
397 scons warp-903334 warp-903336 warp-903338 warp-903342 warp-903344 warp-903346
398 assembleCoadd.py --legacyCoadd $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-R \
399 --selectId visit=903334 ccd=16 --selectId visit=903334 ccd=22 --selectId visit=903334 ccd=23 \
400 --selectId visit=903334 ccd=100 --selectId visit=903336 ccd=17 --selectId visit=903336 ccd=24 \
401 --selectId visit=903338 ccd=18 --selectId visit=903338 ccd=25 --selectId visit=903342 ccd=4 \
402 --selectId visit=903342 ccd=10 --selectId visit=903342 ccd=100 --selectId visit=903344 ccd=0 \
403 --selectId visit=903344 ccd=5 --selectId visit=903344 ccd=11 --selectId visit=903346 ccd=1 \
404 --selectId visit=903346 ccd=6 --selectId visit=903346 ccd=12
406 to generate the coadd for the HSC-R band if you are interested in
407 following multiBand Coadd processing as discussed in `pipeTasks_multiBand`
408 (but note that normally, one would use the `SafeClipAssembleCoaddTask`
409 rather than `AssembleCoaddTask` to make the coadd.
411 ConfigClass = AssembleCoaddConfig
412 _DefaultName =
"assembleCoadd"
414 def __init__(self, *args, **kwargs):
417 argNames = [
"config",
"name",
"parentTask",
"log"]
418 kwargs.update({k: v
for k, v
in zip(argNames, args)})
419 warnings.warn(
"AssembleCoadd received positional args, and casting them as kwargs: %s. "
420 "PipelineTask will not take positional args" % argNames, FutureWarning)
422 super().__init__(**kwargs)
423 self.makeSubtask(
"interpImage")
424 self.makeSubtask(
"scaleZeroPoint")
426 if self.config.doMaskBrightObjects:
429 self.brightObjectBitmask = 1 << mask.addMaskPlane(self.config.brightObjectMaskName)
430 except pexExceptions.LsstCppException:
431 raise RuntimeError(
"Unable to define mask plane for bright objects; planes used are %s" %
432 mask.getMaskPlaneDict().
keys())
435 if self.config.doInputMap:
436 self.makeSubtask(
"inputMapper")
438 self.warpType = self.config.warpType
440 @utils.inheritDoc(pipeBase.PipelineTask)
441 def runQuantum(self, butlerQC, inputRefs, outputRefs):
446 Assemble a coadd from a set of Warps.
448 PipelineTask (Gen3) entry point to Coadd a set of Warps.
449 Analogous to `runDataRef`, it prepares all the data products to be
450 passed to `run`, and processes the results before returning a struct
451 of results to be written out. AssembleCoadd cannot fit all Warps in memory.
452 Therefore, its inputs are accessed subregion by subregion
453 by the Gen3 `DeferredDatasetHandle` that is analagous to the Gen2
454 `lsst.daf.persistence.ButlerDataRef`. Any updates to this method should
455 correspond to an update in `runDataRef` while both entry points
458 inputData = butlerQC.get(inputRefs)
462 skyMap = inputData[
"skyMap"]
463 outputDataId = butlerQC.quantum.dataId
466 tractId=outputDataId[
'tract'],
467 patchId=outputDataId[
'patch'])
469 if self.config.doSelectVisits:
470 warpRefList = self.filterWarps(inputData[
'inputWarps'], inputData[
'selectedVisits'])
472 warpRefList = inputData[
'inputWarps']
475 inputs = self.prepareInputs(warpRefList)
476 self.log.
info(
"Found %d %s", len(inputs.tempExpRefList),
477 self.getTempExpDatasetName(self.warpType))
478 if len(inputs.tempExpRefList) == 0:
479 raise pipeBase.NoWorkFound(
"No coadd temporary exposures found")
481 supplementaryData = self.makeSupplementaryDataGen3(butlerQC, inputRefs, outputRefs)
482 retStruct = self.run(inputData[
'skyInfo'], inputs.tempExpRefList, inputs.imageScalerList,
483 inputs.weightList, supplementaryData=supplementaryData)
485 inputData.setdefault(
'brightObjectMask',
None)
486 self.processResults(retStruct.coaddExposure, inputData[
'brightObjectMask'], outputDataId)
488 if self.config.doWrite:
489 butlerQC.put(retStruct, outputRefs)
493 def runDataRef(self, dataRef, selectDataList=None, warpRefList=None):
494 """Assemble a coadd from a set of Warps.
496 Pipebase.CmdlineTask entry point to Coadd a set of Warps.
497 Compute weights to be applied to each Warp and
498 find scalings to match the photometric zeropoint to a reference Warp.
499 Assemble the Warps using `run`. Interpolate over NaNs and
500 optionally write the coadd to disk. Return the coadded exposure.
504 dataRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
505 Data reference defining the patch for coaddition and the
506 reference Warp (if ``config.autoReference=False``).
507 Used to access the following data products:
508 - ``self.config.coaddName + "Coadd_skyMap"``
509 - ``self.config.coaddName + "Coadd_ + <warpType> + "Warp"`` (optionally)
510 - ``self.config.coaddName + "Coadd"``
511 selectDataList : `list`
512 List of data references to Calexps. Data to be coadded will be
513 selected from this list based on overlap with the patch defined
514 by dataRef, grouped by visit, and converted to a list of data
517 List of data references to Warps to be coadded.
518 Note: `warpRefList` is just the new name for `tempExpRefList`.
522 retStruct : `lsst.pipe.base.Struct`
523 Result struct with components:
525 - ``coaddExposure``: coadded exposure (``Exposure``).
526 - ``nImage``: exposure count image (``Image``).
528 if selectDataList
and warpRefList:
529 raise RuntimeError(
"runDataRef received both a selectDataList and warpRefList, "
530 "and which to use is ambiguous. Please pass only one.")
532 skyInfo = self.getSkyInfo(dataRef)
533 if warpRefList
is None:
534 calExpRefList = self.selectExposures(dataRef, skyInfo, selectDataList=selectDataList)
535 if len(calExpRefList) == 0:
536 self.log.
warning(
"No exposures to coadd")
538 self.log.
info(
"Coadding %d exposures", len(calExpRefList))
540 warpRefList = self.getTempExpRefList(dataRef, calExpRefList)
542 inputData = self.prepareInputs(warpRefList)
543 self.log.
info(
"Found %d %s", len(inputData.tempExpRefList),
544 self.getTempExpDatasetName(self.warpType))
545 if len(inputData.tempExpRefList) == 0:
546 self.log.
warning(
"No coadd temporary exposures found")
549 supplementaryData = self.makeSupplementaryData(dataRef, warpRefList=inputData.tempExpRefList)
551 retStruct = self.run(skyInfo, inputData.tempExpRefList, inputData.imageScalerList,
552 inputData.weightList, supplementaryData=supplementaryData)
554 brightObjects = self.readBrightObjectMasks(dataRef)
if self.config.doMaskBrightObjects
else None
555 self.processResults(retStruct.coaddExposure, brightObjectMasks=brightObjects, dataId=dataRef.dataId)
557 if self.config.doWrite:
558 if self.getCoaddDatasetName(self.warpType) ==
"deepCoadd" and self.config.hasFakes:
559 coaddDatasetName =
"fakes_" + self.getCoaddDatasetName(self.warpType)
561 coaddDatasetName = self.getCoaddDatasetName(self.warpType)
562 self.log.
info(
"Persisting %s", coaddDatasetName)
563 dataRef.put(retStruct.coaddExposure, coaddDatasetName)
564 if self.config.doNImage
and retStruct.nImage
is not None:
565 dataRef.put(retStruct.nImage, self.getCoaddDatasetName(self.warpType) +
'_nImage')
570 """Interpolate over missing data and mask bright stars.
574 coaddExposure : `lsst.afw.image.Exposure`
575 The coadded exposure to process.
576 dataRef : `lsst.daf.persistence.ButlerDataRef`
577 Butler data reference for supplementary data.
579 if self.config.doInterp:
580 self.interpImage.
run(coaddExposure.getMaskedImage(), planeName=
"NO_DATA")
582 varArray = coaddExposure.variance.array
583 with numpy.errstate(invalid=
"ignore"):
584 varArray[:] = numpy.where(varArray > 0, varArray, numpy.inf)
586 if self.config.doMaskBrightObjects:
587 self.setBrightObjectMasks(coaddExposure, brightObjectMasks, dataId)
590 """Make additional inputs to run() specific to subclasses (Gen2)
592 Duplicates interface of `runDataRef` method
593 Available to be implemented by subclasses only if they need the
594 coadd dataRef for performing preliminary processing before
595 assembling the coadd.
599 dataRef : `lsst.daf.persistence.ButlerDataRef`
600 Butler data reference for supplementary data.
601 selectDataList : `list` (optional)
602 Optional List of data references to Calexps.
603 warpRefList : `list` (optional)
604 Optional List of data references to Warps.
606 return pipeBase.Struct()
609 """Make additional inputs to run() specific to subclasses (Gen3)
611 Duplicates interface of `runQuantum` method.
612 Available to be implemented by subclasses only if they need the
613 coadd dataRef for performing preliminary processing before
614 assembling the coadd.
618 butlerQC : `lsst.pipe.base.ButlerQuantumContext`
619 Gen3 Butler object for fetching additional data products before
620 running the Task specialized for quantum being processed
621 inputRefs : `lsst.pipe.base.InputQuantizedConnection`
622 Attributes are the names of the connections describing input dataset types.
623 Values are DatasetRefs that task consumes for corresponding dataset type.
624 DataIds are guaranteed to match data objects in ``inputData``.
625 outputRefs : `lsst.pipe.base.OutputQuantizedConnection`
626 Attributes are the names of the connections describing output dataset types.
627 Values are DatasetRefs that task is to produce
628 for corresponding dataset type.
630 return pipeBase.Struct()
633 """Generate list data references corresponding to warped exposures
634 that lie within the patch to be coadded.
639 Data reference for patch.
640 calExpRefList : `list`
641 List of data references for input calexps.
645 tempExpRefList : `list`
646 List of Warp/CoaddTempExp data references.
648 butler = patchRef.getButler()
649 groupData =
groupPatchExposures(patchRef, calExpRefList, self.getCoaddDatasetName(self.warpType),
650 self.getTempExpDatasetName(self.warpType))
651 tempExpRefList = [
getGroupDataRef(butler, self.getTempExpDatasetName(self.warpType),
652 g, groupData.keys)
for
653 g
in groupData.groups.keys()]
654 return tempExpRefList
657 """Prepare the input warps for coaddition by measuring the weight for
658 each warp and the scaling for the photometric zero point.
660 Each Warp has its own photometric zeropoint and background variance.
661 Before coadding these Warps together, compute a scale factor to
662 normalize the photometric zeropoint and compute the weight for each Warp.
667 List of data references to tempExp
671 result : `lsst.pipe.base.Struct`
672 Result struct with components:
674 - ``tempExprefList``: `list` of data references to tempExp.
675 - ``weightList``: `list` of weightings.
676 - ``imageScalerList``: `list` of image scalers.
679 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
680 statsCtrl.setNumIter(self.config.clipIter)
681 statsCtrl.setAndMask(self.getBadPixelMask())
682 statsCtrl.setNanSafe(
True)
689 tempExpName = self.getTempExpDatasetName(self.warpType)
690 for tempExpRef
in refList:
693 if not isinstance(tempExpRef, DeferredDatasetHandle):
694 if not tempExpRef.datasetExists(tempExpName):
695 self.log.
warning(
"Could not find %s %s; skipping it", tempExpName, tempExpRef.dataId)
698 tempExp = tempExpRef.get(datasetType=tempExpName, immediate=
True)
700 if numpy.isnan(tempExp.image.array).
all():
702 maskedImage = tempExp.getMaskedImage()
703 imageScaler = self.scaleZeroPoint.computeImageScaler(
708 imageScaler.scaleMaskedImage(maskedImage)
709 except Exception
as e:
710 self.log.
warning(
"Scaling failed for %s (skipping it): %s", tempExpRef.dataId, e)
713 afwMath.MEANCLIP, statsCtrl)
714 meanVar, meanVarErr = statObj.getResult(afwMath.MEANCLIP)
715 weight = 1.0 / float(meanVar)
716 if not numpy.isfinite(weight):
717 self.log.
warning(
"Non-finite weight for %s: skipping", tempExpRef.dataId)
719 self.log.
info(
"Weight of %s %s = %0.3f", tempExpName, tempExpRef.dataId, weight)
724 tempExpRefList.append(tempExpRef)
725 weightList.append(weight)
726 imageScalerList.append(imageScaler)
728 return pipeBase.Struct(tempExpRefList=tempExpRefList, weightList=weightList,
729 imageScalerList=imageScalerList)
732 """Prepare the statistics for coadding images.
736 mask : `int`, optional
737 Bit mask value to exclude from coaddition.
741 stats : `lsst.pipe.base.Struct`
742 Statistics structure with the following fields:
744 - ``statsCtrl``: Statistics control object for coadd
745 (`lsst.afw.math.StatisticsControl`)
746 - ``statsFlags``: Statistic for coadd (`lsst.afw.math.Property`)
749 mask = self.getBadPixelMask()
751 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
752 statsCtrl.setNumIter(self.config.clipIter)
753 statsCtrl.setAndMask(mask)
754 statsCtrl.setNanSafe(
True)
755 statsCtrl.setWeighted(
True)
756 statsCtrl.setCalcErrorFromInputVariance(self.config.calcErrorFromInputVariance)
757 for plane, threshold
in self.config.maskPropagationThresholds.items():
758 bit = afwImage.Mask.getMaskPlane(plane)
759 statsCtrl.setMaskPropagationThreshold(bit, threshold)
761 return pipeBase.Struct(ctrl=statsCtrl, flags=statsFlags)
764 def run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
765 altMaskList=None, mask=None, supplementaryData=None):
766 """Assemble a coadd from input warps
768 Assemble the coadd using the provided list of coaddTempExps. Since
769 the full coadd covers a patch (a large area), the assembly is
770 performed over small areas on the image at a time in order to
771 conserve memory usage. Iterate over subregions within the outer
772 bbox of the patch using `assembleSubregion` to stack the corresponding
773 subregions from the coaddTempExps with the statistic specified.
774 Set the edge bits the coadd mask based on the weight map.
778 skyInfo : `lsst.pipe.base.Struct`
779 Struct with geometric information about the patch.
780 tempExpRefList : `list`
781 List of data references to Warps (previously called CoaddTempExps).
782 imageScalerList : `list`
783 List of image scalers.
786 altMaskList : `list`, optional
787 List of alternate masks to use rather than those stored with
789 mask : `int`, optional
790 Bit mask value to exclude from coaddition.
791 supplementaryData : lsst.pipe.base.Struct, optional
792 Struct with additional data products needed to assemble coadd.
793 Only used by subclasses that implement `makeSupplementaryData`
798 result : `lsst.pipe.base.Struct`
799 Result struct with components:
801 - ``coaddExposure``: coadded exposure (``lsst.afw.image.Exposure``).
802 - ``nImage``: exposure count image (``lsst.afw.image.Image``), if requested.
803 - ``inputMap``: bit-wise map of inputs, if requested.
804 - ``warpRefList``: input list of refs to the warps (
805 ``lsst.daf.butler.DeferredDatasetHandle`` or
806 ``lsst.daf.persistence.ButlerDataRef``)
808 - ``imageScalerList``: input list of image scalers (unmodified)
809 - ``weightList``: input list of weights (unmodified)
811 tempExpName = self.getTempExpDatasetName(self.warpType)
812 self.log.
info(
"Assembling %s %s", len(tempExpRefList), tempExpName)
813 stats = self.prepareStats(mask=mask)
815 if altMaskList
is None:
816 altMaskList = [
None]*len(tempExpRefList)
818 coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs)
819 coaddExposure.setPhotoCalib(self.scaleZeroPoint.getPhotoCalib())
820 coaddExposure.getInfo().setCoaddInputs(self.inputRecorder.makeCoaddInputs())
821 self.assembleMetadata(coaddExposure, tempExpRefList, weightList)
822 coaddMaskedImage = coaddExposure.getMaskedImage()
823 subregionSizeArr = self.config.subregionSize
824 subregionSize =
geom.Extent2I(subregionSizeArr[0], subregionSizeArr[1])
826 if self.config.doNImage:
827 nImage = afwImage.ImageU(skyInfo.bbox)
832 if self.config.doInputMap:
833 self.inputMapper.build_ccd_input_map(skyInfo.bbox,
835 coaddExposure.getInfo().getCoaddInputs().ccds)
837 if self.config.doOnlineForMean
and self.config.statistic ==
"MEAN":
839 self.assembleOnlineMeanCoadd(coaddExposure, tempExpRefList, imageScalerList,
840 weightList, altMaskList, stats.ctrl,
842 except Exception
as e:
843 self.log.
fatal(
"Cannot compute online coadd %s", e)
845 for subBBox
in self._subBBoxIter(skyInfo.bbox, subregionSize):
847 self.assembleSubregion(coaddExposure, subBBox, tempExpRefList, imageScalerList,
848 weightList, altMaskList, stats.flags, stats.ctrl,
850 except Exception
as e:
851 self.log.
fatal(
"Cannot compute coadd %s: %s", subBBox, e)
854 if self.config.doInputMap:
855 self.inputMapper.finalize_ccd_input_map_mask()
856 inputMap = self.inputMapper.ccd_input_map
860 self.setInexactPsf(coaddMaskedImage.getMask())
864 return pipeBase.Struct(coaddExposure=coaddExposure, nImage=nImage,
865 warpRefList=tempExpRefList, imageScalerList=imageScalerList,
866 weightList=weightList, inputMap=inputMap)
869 """Set the metadata for the coadd.
871 This basic implementation sets the filter from the first input.
875 coaddExposure : `lsst.afw.image.Exposure`
876 The target exposure for the coadd.
877 tempExpRefList : `list`
878 List of data references to tempExp.
882 assert len(tempExpRefList) == len(weightList),
"Length mismatch"
883 tempExpName = self.getTempExpDatasetName(self.warpType)
889 if isinstance(tempExpRefList[0], DeferredDatasetHandle):
891 tempExpList = [tempExpRef.get(parameters={
'bbox': bbox})
for tempExpRef
in tempExpRefList]
894 tempExpList = [tempExpRef.get(tempExpName +
"_sub", bbox=bbox, immediate=
True)
895 for tempExpRef
in tempExpRefList]
896 numCcds = sum(len(tempExp.getInfo().getCoaddInputs().ccds)
for tempExp
in tempExpList)
901 coaddInputs = coaddExposure.getInfo().getCoaddInputs()
902 coaddInputs.ccds.reserve(numCcds)
903 coaddInputs.visits.reserve(len(tempExpList))
905 for tempExp, weight
in zip(tempExpList, weightList):
906 self.inputRecorder.addVisitToCoadd(coaddInputs, tempExp, weight)
908 if self.config.doUsePsfMatchedPolygons:
909 self.shrinkValidPolygons(coaddInputs)
911 coaddInputs.visits.sort()
912 if self.warpType ==
"psfMatched":
917 modelPsfList = [tempExp.getPsf()
for tempExp
in tempExpList]
918 modelPsfWidthList = [modelPsf.computeBBox().getWidth()
for modelPsf
in modelPsfList]
919 psf = modelPsfList[modelPsfWidthList.index(
max(modelPsfWidthList))]
921 psf = measAlg.CoaddPsf(coaddInputs.ccds, coaddExposure.getWcs(),
922 self.config.coaddPsf.makeControl())
923 coaddExposure.setPsf(psf)
924 apCorrMap = measAlg.makeCoaddApCorrMap(coaddInputs.ccds, coaddExposure.getBBox(afwImage.PARENT),
925 coaddExposure.getWcs())
926 coaddExposure.getInfo().setApCorrMap(apCorrMap)
927 if self.config.doAttachTransmissionCurve:
928 transmissionCurve = measAlg.makeCoaddTransmissionCurve(coaddExposure.getWcs(), coaddInputs.ccds)
929 coaddExposure.getInfo().setTransmissionCurve(transmissionCurve)
932 altMaskList, statsFlags, statsCtrl, nImage=None):
933 """Assemble the coadd for a sub-region.
935 For each coaddTempExp, check for (and swap in) an alternative mask
936 if one is passed. Remove mask planes listed in
937 `config.removeMaskPlanes`. Finally, stack the actual exposures using
938 `lsst.afw.math.statisticsStack` with the statistic specified by
939 statsFlags. Typically, the statsFlag will be one of lsst.afw.math.MEAN for
940 a mean-stack or `lsst.afw.math.MEANCLIP` for outlier rejection using
941 an N-sigma clipped mean where N and iterations are specified by
942 statsCtrl. Assign the stacked subregion back to the coadd.
946 coaddExposure : `lsst.afw.image.Exposure`
947 The target exposure for the coadd.
948 bbox : `lsst.geom.Box`
950 tempExpRefList : `list`
951 List of data reference to tempExp.
952 imageScalerList : `list`
953 List of image scalers.
957 List of alternate masks to use rather than those stored with
958 tempExp, or None. Each element is dict with keys = mask plane
959 name to which to add the spans.
960 statsFlags : `lsst.afw.math.Property`
961 Property object for statistic for coadd.
962 statsCtrl : `lsst.afw.math.StatisticsControl`
963 Statistics control object for coadd.
964 nImage : `lsst.afw.image.ImageU`, optional
965 Keeps track of exposure count for each pixel.
967 self.log.
debug(
"Computing coadd over %s", bbox)
968 tempExpName = self.getTempExpDatasetName(self.warpType)
969 coaddExposure.mask.addMaskPlane(
"REJECTED")
970 coaddExposure.mask.addMaskPlane(
"CLIPPED")
971 coaddExposure.mask.addMaskPlane(
"SENSOR_EDGE")
972 maskMap = self.setRejectedMaskMapping(statsCtrl)
973 clipped = afwImage.Mask.getPlaneBitMask(
"CLIPPED")
975 if nImage
is not None:
976 subNImage = afwImage.ImageU(bbox.getWidth(), bbox.getHeight())
977 for tempExpRef, imageScaler, altMask
in zip(tempExpRefList, imageScalerList, altMaskList):
979 if isinstance(tempExpRef, DeferredDatasetHandle):
981 exposure = tempExpRef.get(parameters={
'bbox': bbox})
984 exposure = tempExpRef.get(tempExpName +
"_sub", bbox=bbox)
986 maskedImage = exposure.getMaskedImage()
987 mask = maskedImage.getMask()
988 if altMask
is not None:
989 self.applyAltMaskPlanes(mask, altMask)
990 imageScaler.scaleMaskedImage(maskedImage)
994 if nImage
is not None:
995 subNImage.getArray()[maskedImage.getMask().getArray() & statsCtrl.getAndMask() == 0] += 1
996 if self.config.removeMaskPlanes:
997 self.removeMaskPlanes(maskedImage)
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"):
1008 coaddExposure.maskedImage.assign(coaddSubregion, bbox)
1009 if nImage
is not None:
1010 nImage.assign(subNImage, bbox)
1013 altMaskList, statsCtrl, nImage=None):
1014 """Assemble the coadd using the "online" method.
1016 This method takes a running sum of images and weights to save memory.
1017 It only works for MEAN statistics.
1021 coaddExposure : `lsst.afw.image.Exposure`
1022 The target exposure for the coadd.
1023 tempExpRefList : `list`
1024 List of data reference to tempExp.
1025 imageScalerList : `list`
1026 List of image scalers.
1029 altMaskList : `list`
1030 List of alternate masks to use rather than those stored with
1031 tempExp, or None. Each element is dict with keys = mask plane
1032 name to which to add the spans.
1033 statsCtrl : `lsst.afw.math.StatisticsControl`
1034 Statistics control object for coadd
1035 nImage : `lsst.afw.image.ImageU`, optional
1036 Keeps track of exposure count for each pixel.
1038 self.log.
debug(
"Computing online coadd.")
1039 tempExpName = self.getTempExpDatasetName(self.warpType)
1040 coaddExposure.mask.addMaskPlane(
"REJECTED")
1041 coaddExposure.mask.addMaskPlane(
"CLIPPED")
1042 coaddExposure.mask.addMaskPlane(
"SENSOR_EDGE")
1043 maskMap = self.setRejectedMaskMapping(statsCtrl)
1044 thresholdDict = AccumulatorMeanStack.stats_ctrl_to_threshold_dict(statsCtrl)
1046 bbox = coaddExposure.maskedImage.getBBox()
1049 coaddExposure.image.array.shape,
1050 statsCtrl.getAndMask(),
1051 mask_threshold_dict=thresholdDict,
1053 no_good_pixels_mask=statsCtrl.getNoGoodPixelsMask(),
1054 calc_error_from_input_variance=self.config.calcErrorFromInputVariance,
1055 compute_n_image=(nImage
is not None)
1058 for tempExpRef, imageScaler, altMask, weight
in zip(tempExpRefList,
1062 if isinstance(tempExpRef, DeferredDatasetHandle):
1064 exposure = tempExpRef.get()
1067 exposure = tempExpRef.get(tempExpName)
1069 maskedImage = exposure.getMaskedImage()
1070 mask = maskedImage.getMask()
1071 if altMask
is not None:
1072 self.applyAltMaskPlanes(mask, altMask)
1073 imageScaler.scaleMaskedImage(maskedImage)
1074 if self.config.removeMaskPlanes:
1075 self.removeMaskPlanes(maskedImage)
1077 stacker.add_masked_image(maskedImage, weight=weight)
1079 if self.config.doInputMap:
1080 visit = exposure.getInfo().getCoaddInputs().visits[0].getId()
1081 self.inputMapper.mask_warp_bbox(bbox, visit, mask, statsCtrl.getAndMask())
1083 stacker.fill_stacked_masked_image(coaddExposure.maskedImage)
1085 if nImage
is not None:
1086 nImage.array[:, :] = stacker.n_image
1089 """Unset the mask of an image for mask planes specified in the config.
1093 maskedImage : `lsst.afw.image.MaskedImage`
1094 The masked image to be modified.
1096 mask = maskedImage.getMask()
1097 for maskPlane
in self.config.removeMaskPlanes:
1099 mask &= ~mask.getPlaneBitMask(maskPlane)
1101 self.log.
debug(
"Unable to remove mask plane %s: no mask plane with that name was found.",
1105 def setRejectedMaskMapping(statsCtrl):
1106 """Map certain mask planes of the warps to new planes for the coadd.
1108 If a pixel is rejected due to a mask value other than EDGE, NO_DATA,
1109 or CLIPPED, set it to REJECTED on the coadd.
1110 If a pixel is rejected due to EDGE, set the coadd pixel to SENSOR_EDGE.
1111 If a pixel is rejected due to CLIPPED, set the coadd pixel to CLIPPED.
1115 statsCtrl : `lsst.afw.math.StatisticsControl`
1116 Statistics control object for coadd
1120 maskMap : `list` of `tuple` of `int`
1121 A list of mappings of mask planes of the warped exposures to
1122 mask planes of the coadd.
1124 edge = afwImage.Mask.getPlaneBitMask(
"EDGE")
1125 noData = afwImage.Mask.getPlaneBitMask(
"NO_DATA")
1126 clipped = afwImage.Mask.getPlaneBitMask(
"CLIPPED")
1127 toReject = statsCtrl.getAndMask() & (~noData) & (~edge) & (~clipped)
1128 maskMap = [(toReject, afwImage.Mask.getPlaneBitMask(
"REJECTED")),
1129 (edge, afwImage.Mask.getPlaneBitMask(
"SENSOR_EDGE")),
1134 """Apply in place alt mask formatted as SpanSets to a mask.
1138 mask : `lsst.afw.image.Mask`
1140 altMaskSpans : `dict`
1141 SpanSet lists to apply. Each element contains the new mask
1142 plane name (e.g. "CLIPPED and/or "NO_DATA") as the key,
1143 and list of SpanSets to apply to the mask.
1147 mask : `lsst.afw.image.Mask`
1150 if self.config.doUsePsfMatchedPolygons:
1151 if (
"NO_DATA" in altMaskSpans)
and (
"NO_DATA" in self.config.badMaskPlanes):
1156 for spanSet
in altMaskSpans[
'NO_DATA']:
1157 spanSet.clippedTo(mask.getBBox()).clearMask(mask, self.getBadPixelMask())
1159 for plane, spanSetList
in altMaskSpans.items():
1160 maskClipValue = mask.addMaskPlane(plane)
1161 for spanSet
in spanSetList:
1162 spanSet.clippedTo(mask.getBBox()).setMask(mask, 2**maskClipValue)
1166 """Shrink coaddInputs' ccds' ValidPolygons in place.
1168 Either modify each ccd's validPolygon in place, or if CoaddInputs
1169 does not have a validPolygon, create one from its bbox.
1173 coaddInputs : `lsst.afw.image.coaddInputs`
1177 for ccd
in coaddInputs.ccds:
1178 polyOrig = ccd.getValidPolygon()
1179 validPolyBBox = polyOrig.getBBox()
if polyOrig
else ccd.getBBox()
1180 validPolyBBox.grow(-self.config.matchingKernelSize//2)
1182 validPolygon = polyOrig.intersectionSingle(validPolyBBox)
1184 validPolygon = afwGeom.polygon.Polygon(
geom.Box2D(validPolyBBox))
1185 ccd.setValidPolygon(validPolygon)
1188 """Retrieve the bright object masks.
1190 Returns None on failure.
1194 dataRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
1199 result : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
1200 Bright object mask from the Butler object, or None if it cannot
1204 return dataRef.get(datasetType=
"brightObjectMask", immediate=
True)
1205 except Exception
as e:
1206 self.log.
warning(
"Unable to read brightObjectMask for %s: %s", dataRef.dataId, e)
1210 """Set the bright object masks.
1214 exposure : `lsst.afw.image.Exposure`
1215 Exposure under consideration.
1216 dataId : `lsst.daf.persistence.dataId`
1217 Data identifier dict for patch.
1218 brightObjectMasks : `lsst.afw.table`
1219 Table of bright objects to mask.
1222 if brightObjectMasks
is None:
1223 self.log.
warning(
"Unable to apply bright object mask: none supplied")
1225 self.log.
info(
"Applying %d bright object masks to %s", len(brightObjectMasks), dataId)
1226 mask = exposure.getMaskedImage().getMask()
1227 wcs = exposure.getWcs()
1228 plateScale = wcs.getPixelScale().asArcseconds()
1230 for rec
in brightObjectMasks:
1231 center =
geom.PointI(wcs.skyToPixel(rec.getCoord()))
1232 if rec[
"type"] ==
"box":
1233 assert rec[
"angle"] == 0.0, (
"Angle != 0 for mask object %s" % rec[
"id"])
1234 width = rec[
"width"].asArcseconds()/plateScale
1235 height = rec[
"height"].asArcseconds()/plateScale
1238 bbox =
geom.Box2I(center - halfSize, center + halfSize)
1241 geom.PointI(int(center[0] + 0.5*width), int(center[1] + 0.5*height)))
1243 elif rec[
"type"] ==
"circle":
1244 radius = int(rec[
"radius"].asArcseconds()/plateScale)
1245 spans = afwGeom.SpanSet.fromShape(radius, offset=center)
1247 self.log.
warning(
"Unexpected region type %s at %s", rec[
"type"], center)
1249 spans.clippedTo(mask.getBBox()).setMask(mask, self.brightObjectBitmask)
1252 """Set INEXACT_PSF mask plane.
1254 If any of the input images isn't represented in the coadd (due to
1255 clipped pixels or chip gaps), the `CoaddPsf` will be inexact. Flag
1260 mask : `lsst.afw.image.Mask`
1261 Coadded exposure's mask, modified in-place.
1263 mask.addMaskPlane(
"INEXACT_PSF")
1264 inexactPsf = mask.getPlaneBitMask(
"INEXACT_PSF")
1265 sensorEdge = mask.getPlaneBitMask(
"SENSOR_EDGE")
1266 clipped = mask.getPlaneBitMask(
"CLIPPED")
1267 rejected = mask.getPlaneBitMask(
"REJECTED")
1268 array = mask.getArray()
1269 selected = array & (sensorEdge | clipped | rejected) > 0
1270 array[selected] |= inexactPsf
1273 def _makeArgumentParser(cls):
1274 """Create an argument parser.
1276 parser = pipeBase.ArgumentParser(name=cls._DefaultName)
1277 parser.add_id_argument(
"--id", cls.ConfigClass().coaddName +
"Coadd_"
1278 + cls.ConfigClass().warpType +
"Warp",
1279 help=
"data ID, e.g. --id tract=12345 patch=1,2",
1280 ContainerClass=AssembleCoaddDataIdContainer)
1281 parser.add_id_argument(
"--selectId",
"calexp", help=
"data ID, e.g. --selectId visit=6789 ccd=0..9",
1282 ContainerClass=SelectDataIdContainer)
1286 def _subBBoxIter(bbox, subregionSize):
1287 """Iterate over subregions of a bbox.
1291 bbox : `lsst.geom.Box2I`
1292 Bounding box over which to iterate.
1293 subregionSize: `lsst.geom.Extent2I`
1298 subBBox : `lsst.geom.Box2I`
1299 Next sub-bounding box of size ``subregionSize`` or smaller; each ``subBBox``
1300 is contained within ``bbox``, so it may be smaller than ``subregionSize`` at
1301 the edges of ``bbox``, but it will never be empty.
1304 raise RuntimeError(
"bbox %s is empty" % (bbox,))
1305 if subregionSize[0] < 1
or subregionSize[1] < 1:
1306 raise RuntimeError(
"subregionSize %s must be nonzero" % (subregionSize,))
1308 for rowShift
in range(0, bbox.getHeight(), subregionSize[1]):
1309 for colShift
in range(0, bbox.getWidth(), subregionSize[0]):
1312 if subBBox.isEmpty():
1313 raise RuntimeError(
"Bug: empty bbox! bbox=%s, subregionSize=%s, "
1314 "colShift=%s, rowShift=%s" %
1315 (bbox, subregionSize, colShift, rowShift))
1319 """Return list of only inputRefs with visitId in goodVisits ordered by goodVisit
1324 List of `lsst.pipe.base.connections.DeferredDatasetRef` with dataId containing visit
1326 Dictionary with good visitIds as the keys. Value ignored.
1330 filteredInputs : `list`
1331 Filtered and sorted list of `lsst.pipe.base.connections.DeferredDatasetRef`
1333 inputWarpDict = {inputRef.ref.dataId[
'visit']: inputRef
for inputRef
in inputs}
1335 for visit
in goodVisits.keys():
1336 if visit
in inputWarpDict:
1337 filteredInputs.append(inputWarpDict[visit])
1338 return filteredInputs
1342 """A version of `lsst.pipe.base.DataIdContainer` specialized for assembleCoadd.
1346 """Make self.refList from self.idList.
1351 Results of parsing command-line (with ``butler`` and ``log`` elements).
1353 datasetType = namespace.config.coaddName +
"Coadd"
1354 keysCoadd = namespace.butler.getKeys(datasetType=datasetType, level=self.level)
1356 for dataId
in self.idList:
1358 for key
in keysCoadd:
1359 if key
not in dataId:
1360 raise RuntimeError(
"--id must include " + key)
1362 dataRef = namespace.butler.dataRef(
1363 datasetType=datasetType,
1366 self.refList.
append(dataRef)
1370 """Function to count the number of pixels with a specific mask in a
1373 Find the intersection of mask & footprint. Count all pixels in the mask
1374 that are in the intersection that have bitmask set but do not have
1375 ignoreMask set. Return the count.
1379 mask : `lsst.afw.image.Mask`
1380 Mask to define intersection region by.
1381 footprint : `lsst.afw.detection.Footprint`
1382 Footprint to define the intersection region by.
1384 Specific mask that we wish to count the number of occurances of.
1386 Pixels to not consider.
1391 Count of number of pixels in footprint with specified mask.
1393 bbox = footprint.getBBox()
1394 bbox.clip(mask.getBBox(afwImage.PARENT))
1396 subMask = mask.Factory(mask, bbox, afwImage.PARENT)
1397 footprint.spans.setMask(fp, bitmask)
1398 return numpy.logical_and((subMask.getArray() & fp.getArray()) > 0,
1399 (subMask.getArray() & ignoreMask) == 0).sum()
1403 """Configuration parameters for the SafeClipAssembleCoaddTask.
1405 clipDetection = pexConfig.ConfigurableField(
1406 target=SourceDetectionTask,
1407 doc=
"Detect sources on difference between unclipped and clipped coadd")
1408 minClipFootOverlap = pexConfig.Field(
1409 doc=
"Minimum fractional overlap of clipped footprint with visit DETECTED to be clipped",
1413 minClipFootOverlapSingle = pexConfig.Field(
1414 doc=
"Minimum fractional overlap of clipped footprint with visit DETECTED to be "
1415 "clipped when only one visit overlaps",
1419 minClipFootOverlapDouble = pexConfig.Field(
1420 doc=
"Minimum fractional overlap of clipped footprints with visit DETECTED to be "
1421 "clipped when two visits overlap",
1425 maxClipFootOverlapDouble = pexConfig.Field(
1426 doc=
"Maximum fractional overlap of clipped footprints with visit DETECTED when "
1427 "considering two visits",
1431 minBigOverlap = pexConfig.Field(
1432 doc=
"Minimum number of pixels in footprint to use DETECTED mask from the single visits "
1433 "when labeling clipped footprints",
1439 """Set default values for clipDetection.
1443 The numeric values for these configuration parameters were
1444 empirically determined, future work may further refine them.
1446 AssembleCoaddConfig.setDefaults(self)
1447 self.
clipDetectionclipDetection.doTempLocalBackground =
False
1448 self.
clipDetectionclipDetection.reEstimateBackground =
False
1449 self.
clipDetectionclipDetection.returnOriginalFootprints =
False
1455 self.
clipDetectionclipDetection.thresholdType =
"pixel_stdev"
1462 log.warning(
"Additional Sigma-clipping not allowed in Safe-clipped Coadds. "
1463 "Ignoring doSigmaClip.")
1466 raise ValueError(
"Only MEAN statistic allowed for final stacking in SafeClipAssembleCoadd "
1467 "(%s chosen). Please set statistic to MEAN."
1469 AssembleCoaddTask.ConfigClass.validate(self)
1473 """Assemble a coadded image from a set of coadded temporary exposures,
1474 being careful to clip & flag areas with potential artifacts.
1476 In ``AssembleCoaddTask``, we compute the coadd as an clipped mean (i.e.,
1477 we clip outliers). The problem with doing this is that when computing the
1478 coadd PSF at a given location, individual visit PSFs from visits with
1479 outlier pixels contribute to the coadd PSF and cannot be treated correctly.
1480 In this task, we correct for this behavior by creating a new
1481 ``badMaskPlane`` 'CLIPPED'. We populate this plane on the input
1482 coaddTempExps and the final coadd where
1484 i. difference imaging suggests that there is an outlier and
1485 ii. this outlier appears on only one or two images.
1487 Such regions will not contribute to the final coadd. Furthermore, any
1488 routine to determine the coadd PSF can now be cognizant of clipped regions.
1489 Note that the algorithm implemented by this task is preliminary and works
1490 correctly for HSC data. Parameter modifications and or considerable
1491 redesigning of the algorithm is likley required for other surveys.
1493 ``SafeClipAssembleCoaddTask`` uses a ``SourceDetectionTask``
1494 "clipDetection" subtask and also sub-classes ``AssembleCoaddTask``.
1495 You can retarget the ``SourceDetectionTask`` "clipDetection" subtask
1500 The `lsst.pipe.base.cmdLineTask.CmdLineTask` interface supports a
1501 flag ``-d`` to import ``debug.py`` from your ``PYTHONPATH``;
1502 see `baseDebug` for more about ``debug.py`` files.
1503 `SafeClipAssembleCoaddTask` has no debug variables of its own.
1504 The ``SourceDetectionTask`` "clipDetection" subtasks may support debug
1505 variables. See the documetation for `SourceDetectionTask` "clipDetection"
1506 for further information.
1510 `SafeClipAssembleCoaddTask` assembles a set of warped ``coaddTempExp``
1511 images into a coadded image. The `SafeClipAssembleCoaddTask` is invoked by
1512 running assembleCoadd.py *without* the flag '--legacyCoadd'.
1514 Usage of ``assembleCoadd.py`` expects a data reference to the tract patch
1515 and filter to be coadded (specified using
1516 '--id = [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]')
1517 along with a list of coaddTempExps to attempt to coadd (specified using
1518 '--selectId [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]').
1519 Only the coaddTempExps that cover the specified tract and patch will be
1520 coadded. A list of the available optional arguments can be obtained by
1521 calling assembleCoadd.py with the --help command line argument:
1523 .. code-block:: none
1525 assembleCoadd.py --help
1527 To demonstrate usage of the `SafeClipAssembleCoaddTask` in the larger
1528 context of multi-band processing, we will generate the HSC-I & -R band
1529 coadds from HSC engineering test data provided in the ci_hsc package.
1530 To begin, assuming that the lsst stack has been already set up, we must
1531 set up the obs_subaru and ci_hsc packages. This defines the environment
1532 variable $CI_HSC_DIR and points at the location of the package. The raw
1533 HSC data live in the ``$CI_HSC_DIR/raw`` directory. To begin assembling
1534 the coadds, we must first
1537 process the individual ccds in $CI_HSC_RAW to produce calibrated exposures
1539 create a skymap that covers the area of the sky present in the raw exposures
1540 - ``makeCoaddTempExp``
1541 warp the individual calibrated exposures to the tangent plane of the coadd</DD>
1543 We can perform all of these steps by running
1545 .. code-block:: none
1547 $CI_HSC_DIR scons warp-903986 warp-904014 warp-903990 warp-904010 warp-903988
1549 This will produce warped coaddTempExps for each visit. To coadd the
1550 warped data, we call ``assembleCoadd.py`` as follows:
1552 .. code-block:: none
1554 assembleCoadd.py $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-I \
1555 --selectId visit=903986 ccd=16 --selectId visit=903986 ccd=22 --selectId visit=903986 ccd=23 \
1556 --selectId visit=903986 ccd=100--selectId visit=904014 ccd=1 --selectId visit=904014 ccd=6 \
1557 --selectId visit=904014 ccd=12 --selectId visit=903990 ccd=18 --selectId visit=903990 ccd=25 \
1558 --selectId visit=904010 ccd=4 --selectId visit=904010 ccd=10 --selectId visit=904010 ccd=100 \
1559 --selectId visit=903988 ccd=16 --selectId visit=903988 ccd=17 --selectId visit=903988 ccd=23 \
1560 --selectId visit=903988 ccd=24
1562 This will process the HSC-I band data. The results are written in
1563 ``$CI_HSC_DIR/DATA/deepCoadd-results/HSC-I``.
1565 You may also choose to run:
1567 .. code-block:: none
1569 scons warp-903334 warp-903336 warp-903338 warp-903342 warp-903344 warp-903346 nnn
1570 assembleCoadd.py $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-R --selectId visit=903334 ccd=16 \
1571 --selectId visit=903334 ccd=22 --selectId visit=903334 ccd=23 --selectId visit=903334 ccd=100 \
1572 --selectId visit=903336 ccd=17 --selectId visit=903336 ccd=24 --selectId visit=903338 ccd=18 \
1573 --selectId visit=903338 ccd=25 --selectId visit=903342 ccd=4 --selectId visit=903342 ccd=10 \
1574 --selectId visit=903342 ccd=100 --selectId visit=903344 ccd=0 --selectId visit=903344 ccd=5 \
1575 --selectId visit=903344 ccd=11 --selectId visit=903346 ccd=1 --selectId visit=903346 ccd=6 \
1576 --selectId visit=903346 ccd=12
1578 to generate the coadd for the HSC-R band if you are interested in following
1579 multiBand Coadd processing as discussed in ``pipeTasks_multiBand``.
1581 ConfigClass = SafeClipAssembleCoaddConfig
1582 _DefaultName =
"safeClipAssembleCoadd"
1585 AssembleCoaddTask.__init__(self, *args, **kwargs)
1586 schema = afwTable.SourceTable.makeMinimalSchema()
1587 self.makeSubtask(
"clipDetection", schema=schema)
1589 @utils.inheritDoc(AssembleCoaddTask)
1591 def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, *args, **kwargs):
1592 """Assemble the coadd for a region.
1594 Compute the difference of coadds created with and without outlier
1595 rejection to identify coadd pixels that have outlier values in some
1597 Detect clipped regions on the difference image and mark these regions
1598 on the one or two individual coaddTempExps where they occur if there
1599 is significant overlap between the clipped region and a source. This
1600 leaves us with a set of footprints from the difference image that have
1601 been identified as having occured on just one or two individual visits.
1602 However, these footprints were generated from a difference image. It
1603 is conceivable for a large diffuse source to have become broken up
1604 into multiple footprints acrosss the coadd difference in this process.
1605 Determine the clipped region from all overlapping footprints from the
1606 detected sources in each visit - these are big footprints.
1607 Combine the small and big clipped footprints and mark them on a new
1609 Generate the coadd using `AssembleCoaddTask.run` without outlier
1610 removal. Clipped footprints will no longer make it into the coadd
1611 because they are marked in the new bad mask plane.
1615 args and kwargs are passed but ignored in order to match the call
1616 signature expected by the parent task.
1618 exp = self.
buildDifferenceImagebuildDifferenceImage(skyInfo, tempExpRefList, imageScalerList, weightList)
1619 mask = exp.getMaskedImage().getMask()
1620 mask.addMaskPlane(
"CLIPPED")
1622 result = self.
detectClipdetectClip(exp, tempExpRefList)
1624 self.log.
info(
'Found %d clipped objects', len(result.clipFootprints))
1626 maskClipValue = mask.getPlaneBitMask(
"CLIPPED")
1627 maskDetValue = mask.getPlaneBitMask(
"DETECTED") | mask.getPlaneBitMask(
"DETECTED_NEGATIVE")
1629 bigFootprints = self.
detectClipBigdetectClipBig(result.clipSpans, result.clipFootprints, result.clipIndices,
1630 result.detectionFootprints, maskClipValue, maskDetValue,
1633 maskClip = mask.Factory(mask.getBBox(afwImage.PARENT))
1634 afwDet.setMaskFromFootprintList(maskClip, result.clipFootprints, maskClipValue)
1636 maskClipBig = maskClip.Factory(mask.getBBox(afwImage.PARENT))
1637 afwDet.setMaskFromFootprintList(maskClipBig, bigFootprints, maskClipValue)
1638 maskClip |= maskClipBig
1641 badMaskPlanes = self.config.badMaskPlanes[:]
1642 badMaskPlanes.append(
"CLIPPED")
1643 badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes)
1644 return AssembleCoaddTask.run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
1645 result.clipSpans, mask=badPixelMask)
1648 """Return an exposure that contains the difference between unclipped
1651 Generate a difference image between clipped and unclipped coadds.
1652 Compute the difference image by subtracting an outlier-clipped coadd
1653 from an outlier-unclipped coadd. Return the difference image.
1657 skyInfo : `lsst.pipe.base.Struct`
1658 Patch geometry information, from getSkyInfo
1659 tempExpRefList : `list`
1660 List of data reference to tempExp
1661 imageScalerList : `list`
1662 List of image scalers
1668 exp : `lsst.afw.image.Exposure`
1669 Difference image of unclipped and clipped coadd wrapped in an Exposure
1671 config = AssembleCoaddConfig()
1676 configIntersection = {k: getattr(self.config, k)
1677 for k, v
in self.config.toDict().
items()
1678 if (k
in config.keys()
and k !=
"connections")}
1679 configIntersection[
'doInputMap'] =
False
1680 configIntersection[
'doNImage'] =
False
1681 config.update(**configIntersection)
1684 config.statistic =
'MEAN'
1685 task = AssembleCoaddTask(config=config)
1686 coaddMean = task.run(skyInfo, tempExpRefList, imageScalerList, weightList).coaddExposure
1688 config.statistic =
'MEANCLIP'
1689 task = AssembleCoaddTask(config=config)
1690 coaddClip = task.run(skyInfo, tempExpRefList, imageScalerList, weightList).coaddExposure
1692 coaddDiff = coaddMean.getMaskedImage().
Factory(coaddMean.getMaskedImage())
1693 coaddDiff -= coaddClip.getMaskedImage()
1694 exp = afwImage.ExposureF(coaddDiff)
1695 exp.setPsf(coaddMean.getPsf())
1699 """Detect clipped regions on an exposure and set the mask on the
1700 individual tempExp masks.
1702 Detect footprints in the difference image after smoothing the
1703 difference image with a Gaussian kernal. Identify footprints that
1704 overlap with one or two input ``coaddTempExps`` by comparing the
1705 computed overlap fraction to thresholds set in the config. A different
1706 threshold is applied depending on the number of overlapping visits
1707 (restricted to one or two). If the overlap exceeds the thresholds,
1708 the footprint is considered "CLIPPED" and is marked as such on the
1709 coaddTempExp. Return a struct with the clipped footprints, the indices
1710 of the ``coaddTempExps`` that end up overlapping with the clipped
1711 footprints, and a list of new masks for the ``coaddTempExps``.
1715 exp : `lsst.afw.image.Exposure`
1716 Exposure to run detection on.
1717 tempExpRefList : `list`
1718 List of data reference to tempExp.
1722 result : `lsst.pipe.base.Struct`
1723 Result struct with components:
1725 - ``clipFootprints``: list of clipped footprints.
1726 - ``clipIndices``: indices for each ``clippedFootprint`` in
1728 - ``clipSpans``: List of dictionaries containing spanSet lists
1729 to clip. Each element contains the new maskplane name
1730 ("CLIPPED") as the key and list of ``SpanSets`` as the value.
1731 - ``detectionFootprints``: List of DETECTED/DETECTED_NEGATIVE plane
1732 compressed into footprints.
1734 mask = exp.getMaskedImage().getMask()
1735 maskDetValue = mask.getPlaneBitMask(
"DETECTED") | mask.getPlaneBitMask(
"DETECTED_NEGATIVE")
1736 fpSet = self.clipDetection.detectFootprints(exp, doSmooth=
True, clearMask=
True)
1738 fpSet.positive.merge(fpSet.negative)
1739 footprints = fpSet.positive
1740 self.log.
info(
'Found %d potential clipped objects', len(footprints.getFootprints()))
1741 ignoreMask = self.getBadPixelMask()
1745 artifactSpanSets = [{
'CLIPPED':
list()}
for _
in tempExpRefList]
1748 visitDetectionFootprints = []
1750 dims = [len(tempExpRefList), len(footprints.getFootprints())]
1751 overlapDetArr = numpy.zeros(dims, dtype=numpy.uint16)
1752 ignoreArr = numpy.zeros(dims, dtype=numpy.uint16)
1755 for i, warpRef
in enumerate(tempExpRefList):
1756 tmpExpMask = warpRef.get(datasetType=self.getTempExpDatasetName(self.warpType),
1757 immediate=
True).getMaskedImage().getMask()
1758 maskVisitDet = tmpExpMask.Factory(tmpExpMask, tmpExpMask.getBBox(afwImage.PARENT),
1759 afwImage.PARENT,
True)
1760 maskVisitDet &= maskDetValue
1762 visitDetectionFootprints.append(visitFootprints)
1764 for j, footprint
in enumerate(footprints.getFootprints()):
1769 for j, footprint
in enumerate(footprints.getFootprints()):
1770 nPixel = footprint.getArea()
1773 for i
in range(len(tempExpRefList)):
1774 ignore = ignoreArr[i, j]
1775 overlapDet = overlapDetArr[i, j]
1776 totPixel = nPixel - ignore
1779 if ignore > overlapDet
or totPixel <= 0.5*nPixel
or overlapDet == 0:
1781 overlap.append(overlapDet/float(totPixel))
1784 overlap = numpy.array(overlap)
1785 if not len(overlap):
1792 if len(overlap) == 1:
1793 if overlap[0] > self.config.minClipFootOverlapSingle:
1798 clipIndex = numpy.where(overlap > self.config.minClipFootOverlap)[0]
1799 if len(clipIndex) == 1:
1801 keepIndex = [clipIndex[0]]
1804 clipIndex = numpy.where(overlap > self.config.minClipFootOverlapDouble)[0]
1805 if len(clipIndex) == 2
and len(overlap) > 3:
1806 clipIndexComp = numpy.where(overlap <= self.config.minClipFootOverlapDouble)[0]
1807 if numpy.max(overlap[clipIndexComp]) <= self.config.maxClipFootOverlapDouble:
1809 keepIndex = clipIndex
1814 for index
in keepIndex:
1815 globalIndex = indexList[index]
1816 artifactSpanSets[globalIndex][
'CLIPPED'].
append(footprint.spans)
1818 clipIndices.append(numpy.array(indexList)[keepIndex])
1819 clipFootprints.append(footprint)
1821 return pipeBase.Struct(clipFootprints=clipFootprints, clipIndices=clipIndices,
1822 clipSpans=artifactSpanSets, detectionFootprints=visitDetectionFootprints)
1824 def detectClipBig(self, clipList, clipFootprints, clipIndices, detectionFootprints,
1825 maskClipValue, maskDetValue, coaddBBox):
1826 """Return individual warp footprints for large artifacts and append
1827 them to ``clipList`` in place.
1829 Identify big footprints composed of many sources in the coadd
1830 difference that may have originated in a large diffuse source in the
1831 coadd. We do this by indentifying all clipped footprints that overlap
1832 significantly with each source in all the coaddTempExps.
1837 List of alt mask SpanSets with clipping information. Modified.
1838 clipFootprints : `list`
1839 List of clipped footprints.
1840 clipIndices : `list`
1841 List of which entries in tempExpClipList each footprint belongs to.
1843 Mask value of clipped pixels.
1845 Mask value of detected pixels.
1846 coaddBBox : `lsst.geom.Box`
1847 BBox of the coadd and warps.
1851 bigFootprintsCoadd : `list`
1852 List of big footprints
1854 bigFootprintsCoadd = []
1855 ignoreMask = self.getBadPixelMask()
1856 for index, (clippedSpans, visitFootprints)
in enumerate(zip(clipList, detectionFootprints)):
1857 maskVisitDet = afwImage.MaskX(coaddBBox, 0x0)
1858 for footprint
in visitFootprints.getFootprints():
1859 footprint.spans.setMask(maskVisitDet, maskDetValue)
1862 clippedFootprintsVisit = []
1863 for foot, clipIndex
in zip(clipFootprints, clipIndices):
1864 if index
not in clipIndex:
1866 clippedFootprintsVisit.append(foot)
1867 maskVisitClip = maskVisitDet.Factory(maskVisitDet.getBBox(afwImage.PARENT))
1868 afwDet.setMaskFromFootprintList(maskVisitClip, clippedFootprintsVisit, maskClipValue)
1870 bigFootprintsVisit = []
1871 for foot
in visitFootprints.getFootprints():
1872 if foot.getArea() < self.config.minBigOverlap:
1875 if nCount > self.config.minBigOverlap:
1876 bigFootprintsVisit.append(foot)
1877 bigFootprintsCoadd.append(foot)
1879 for footprint
in bigFootprintsVisit:
1880 clippedSpans[
"CLIPPED"].
append(footprint.spans)
1882 return bigFootprintsCoadd
1886 psfMatchedWarps = pipeBase.connectionTypes.Input(
1887 doc=(
"PSF-Matched Warps are required by CompareWarp regardless of the coadd type requested. "
1888 "Only PSF-Matched Warps make sense for image subtraction. "
1889 "Therefore, they must be an additional declared input."),
1890 name=
"{inputCoaddName}Coadd_psfMatchedWarp",
1891 storageClass=
"ExposureF",
1892 dimensions=(
"tract",
"patch",
"skymap",
"visit"),
1896 templateCoadd = pipeBase.connectionTypes.Output(
1897 doc=(
"Model of the static sky, used to find temporal artifacts. Typically a PSF-Matched, "
1898 "sigma-clipped coadd. Written if and only if assembleStaticSkyModel.doWrite=True"),
1899 name=
"{outputCoaddName}CoaddPsfMatched",
1900 storageClass=
"ExposureF",
1901 dimensions=(
"tract",
"patch",
"skymap",
"band"),
1906 if not config.assembleStaticSkyModel.doWrite:
1907 self.outputs.remove(
"templateCoadd")
1912 pipelineConnections=CompareWarpAssembleCoaddConnections):
1913 assembleStaticSkyModel = pexConfig.ConfigurableField(
1914 target=AssembleCoaddTask,
1915 doc=
"Task to assemble an artifact-free, PSF-matched Coadd to serve as a"
1916 " naive/first-iteration model of the static sky.",
1918 detect = pexConfig.ConfigurableField(
1919 target=SourceDetectionTask,
1920 doc=
"Detect outlier sources on difference between each psfMatched warp and static sky model"
1922 detectTemplate = pexConfig.ConfigurableField(
1923 target=SourceDetectionTask,
1924 doc=
"Detect sources on static sky model. Only used if doPreserveContainedBySource is True"
1926 maskStreaks = pexConfig.ConfigurableField(
1927 target=MaskStreaksTask,
1928 doc=
"Detect streaks on difference between each psfMatched warp and static sky model. Only used if "
1929 "doFilterMorphological is True. Adds a mask plane to an exposure, with the mask plane name set by"
1932 streakMaskName = pexConfig.Field(
1935 doc=
"Name of mask bit used for streaks"
1937 maxNumEpochs = pexConfig.Field(
1938 doc=
"Charactistic maximum local number of epochs/visits in which an artifact candidate can appear "
1939 "and still be masked. The effective maxNumEpochs is a broken linear function of local "
1940 "number of epochs (N): min(maxFractionEpochsLow*N, maxNumEpochs + maxFractionEpochsHigh*N). "
1941 "For each footprint detected on the image difference between the psfMatched warp and static sky "
1942 "model, if a significant fraction of pixels (defined by spatialThreshold) are residuals in more "
1943 "than the computed effective maxNumEpochs, the artifact candidate is deemed persistant rather "
1944 "than transient and not masked.",
1948 maxFractionEpochsLow = pexConfig.RangeField(
1949 doc=
"Fraction of local number of epochs (N) to use as effective maxNumEpochs for low N. "
1950 "Effective maxNumEpochs = "
1951 "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)",
1956 maxFractionEpochsHigh = pexConfig.RangeField(
1957 doc=
"Fraction of local number of epochs (N) to use as effective maxNumEpochs for high N. "
1958 "Effective maxNumEpochs = "
1959 "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)",
1964 spatialThreshold = pexConfig.RangeField(
1965 doc=
"Unitless fraction of pixels defining how much of the outlier region has to meet the "
1966 "temporal criteria. If 0, clip all. If 1, clip none.",
1970 inclusiveMin=
True, inclusiveMax=
True
1972 doScaleWarpVariance = pexConfig.Field(
1973 doc=
"Rescale Warp variance plane using empirical noise?",
1977 scaleWarpVariance = pexConfig.ConfigurableField(
1978 target=ScaleVarianceTask,
1979 doc=
"Rescale variance on warps",
1981 doPreserveContainedBySource = pexConfig.Field(
1982 doc=
"Rescue artifacts from clipping that completely lie within a footprint detected"
1983 "on the PsfMatched Template Coadd. Replicates a behavior of SafeClip.",
1987 doPrefilterArtifacts = pexConfig.Field(
1988 doc=
"Ignore artifact candidates that are mostly covered by the bad pixel mask, "
1989 "because they will be excluded anyway. This prevents them from contributing "
1990 "to the outlier epoch count image and potentially being labeled as persistant."
1991 "'Mostly' is defined by the config 'prefilterArtifactsRatio'.",
1995 prefilterArtifactsMaskPlanes = pexConfig.ListField(
1996 doc=
"Prefilter artifact candidates that are mostly covered by these bad mask planes.",
1998 default=(
'NO_DATA',
'BAD',
'SAT',
'SUSPECT'),
2000 prefilterArtifactsRatio = pexConfig.Field(
2001 doc=
"Prefilter artifact candidates with less than this fraction overlapping good pixels",
2005 doFilterMorphological = pexConfig.Field(
2006 doc=
"Filter artifact candidates based on morphological criteria, i.g. those that appear to "
2013 AssembleCoaddConfig.setDefaults(self)
2019 if "EDGE" in self.badMaskPlanes:
2020 self.badMaskPlanes.remove(
'EDGE')
2021 self.removeMaskPlanes.
append(
'EDGE')
2030 self.
detectdetect.doTempLocalBackground =
False
2031 self.
detectdetect.reEstimateBackground =
False
2032 self.
detectdetect.returnOriginalFootprints =
False
2033 self.
detectdetect.thresholdPolarity =
"both"
2034 self.
detectdetect.thresholdValue = 5
2035 self.
detectdetect.minPixels = 4
2036 self.
detectdetect.isotropicGrow =
True
2037 self.
detectdetect.thresholdType =
"pixel_stdev"
2038 self.
detectdetect.nSigmaToGrow = 0.4
2044 self.
detectTemplatedetectTemplate.returnOriginalFootprints =
False
2049 raise ValueError(
"No dataset type exists for a PSF-Matched Template N Image."
2050 "Please set assembleStaticSkyModel.doNImage=False")
2053 raise ValueError(
"warpType (%s) == assembleStaticSkyModel.warpType (%s) and will compete for "
2054 "the same dataset name. Please set assembleStaticSkyModel.doWrite to False "
2055 "or warpType to 'direct'. assembleStaticSkyModel.warpType should ways be "
2060 """Assemble a compareWarp coadded image from a set of warps
2061 by masking artifacts detected by comparing PSF-matched warps.
2063 In ``AssembleCoaddTask``, we compute the coadd as an clipped mean (i.e.,
2064 we clip outliers). The problem with doing this is that when computing the
2065 coadd PSF at a given location, individual visit PSFs from visits with
2066 outlier pixels contribute to the coadd PSF and cannot be treated correctly.
2067 In this task, we correct for this behavior by creating a new badMaskPlane
2068 'CLIPPED' which marks pixels in the individual warps suspected to contain
2069 an artifact. We populate this plane on the input warps by comparing
2070 PSF-matched warps with a PSF-matched median coadd which serves as a
2071 model of the static sky. Any group of pixels that deviates from the
2072 PSF-matched template coadd by more than config.detect.threshold sigma,
2073 is an artifact candidate. The candidates are then filtered to remove
2074 variable sources and sources that are difficult to subtract such as
2075 bright stars. This filter is configured using the config parameters
2076 ``temporalThreshold`` and ``spatialThreshold``. The temporalThreshold is
2077 the maximum fraction of epochs that the deviation can appear in and still
2078 be considered an artifact. The spatialThreshold is the maximum fraction of
2079 pixels in the footprint of the deviation that appear in other epochs
2080 (where other epochs is defined by the temporalThreshold). If the deviant
2081 region meets this criteria of having a significant percentage of pixels
2082 that deviate in only a few epochs, these pixels have the 'CLIPPED' bit
2083 set in the mask. These regions will not contribute to the final coadd.
2084 Furthermore, any routine to determine the coadd PSF can now be cognizant
2085 of clipped regions. Note that the algorithm implemented by this task is
2086 preliminary and works correctly for HSC data. Parameter modifications and
2087 or considerable redesigning of the algorithm is likley required for other
2090 ``CompareWarpAssembleCoaddTask`` sub-classes
2091 ``AssembleCoaddTask`` and instantiates ``AssembleCoaddTask``
2092 as a subtask to generate the TemplateCoadd (the model of the static sky).
2096 The `lsst.pipe.base.cmdLineTask.CmdLineTask` interface supports a
2097 flag ``-d`` to import ``debug.py`` from your ``PYTHONPATH``; see
2098 ``baseDebug`` for more about ``debug.py`` files.
2100 This task supports the following debug variables:
2103 If True then save the Epoch Count Image as a fits file in the `figPath`
2105 Path to save the debug fits images and figures
2107 For example, put something like:
2109 .. code-block:: python
2112 def DebugInfo(name):
2113 di = lsstDebug.getInfo(name)
2114 if name == "lsst.pipe.tasks.assembleCoadd":
2115 di.saveCountIm = True
2116 di.figPath = "/desired/path/to/debugging/output/images"
2118 lsstDebug.Info = DebugInfo
2120 into your ``debug.py`` file and run ``assemebleCoadd.py`` with the
2121 ``--debug`` flag. Some subtasks may have their own debug variables;
2122 see individual Task documentation.
2126 ``CompareWarpAssembleCoaddTask`` assembles a set of warped images into a
2127 coadded image. The ``CompareWarpAssembleCoaddTask`` is invoked by running
2128 ``assembleCoadd.py`` with the flag ``--compareWarpCoadd``.
2129 Usage of ``assembleCoadd.py`` expects a data reference to the tract patch
2130 and filter to be coadded (specified using
2131 '--id = [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]')
2132 along with a list of coaddTempExps to attempt to coadd (specified using
2133 '--selectId [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]').
2134 Only the warps that cover the specified tract and patch will be coadded.
2135 A list of the available optional arguments can be obtained by calling
2136 ``assembleCoadd.py`` with the ``--help`` command line argument:
2138 .. code-block:: none
2140 assembleCoadd.py --help
2142 To demonstrate usage of the ``CompareWarpAssembleCoaddTask`` in the larger
2143 context of multi-band processing, we will generate the HSC-I & -R band
2144 oadds from HSC engineering test data provided in the ``ci_hsc`` package.
2145 To begin, assuming that the lsst stack has been already set up, we must
2146 set up the ``obs_subaru`` and ``ci_hsc`` packages.
2147 This defines the environment variable ``$CI_HSC_DIR`` and points at the
2148 location of the package. The raw HSC data live in the ``$CI_HSC_DIR/raw``
2149 directory. To begin assembling the coadds, we must first
2152 process the individual ccds in $CI_HSC_RAW to produce calibrated exposures
2154 create a skymap that covers the area of the sky present in the raw exposures
2156 warp the individual calibrated exposures to the tangent plane of the coadd
2158 We can perform all of these steps by running
2160 .. code-block:: none
2162 $CI_HSC_DIR scons warp-903986 warp-904014 warp-903990 warp-904010 warp-903988
2164 This will produce warped ``coaddTempExps`` for each visit. To coadd the
2165 warped data, we call ``assembleCoadd.py`` as follows:
2167 .. code-block:: none
2169 assembleCoadd.py --compareWarpCoadd $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-I \
2170 --selectId visit=903986 ccd=16 --selectId visit=903986 ccd=22 --selectId visit=903986 ccd=23 \
2171 --selectId visit=903986 ccd=100 --selectId visit=904014 ccd=1 --selectId visit=904014 ccd=6 \
2172 --selectId visit=904014 ccd=12 --selectId visit=903990 ccd=18 --selectId visit=903990 ccd=25 \
2173 --selectId visit=904010 ccd=4 --selectId visit=904010 ccd=10 --selectId visit=904010 ccd=100 \
2174 --selectId visit=903988 ccd=16 --selectId visit=903988 ccd=17 --selectId visit=903988 ccd=23 \
2175 --selectId visit=903988 ccd=24
2177 This will process the HSC-I band data. The results are written in
2178 ``$CI_HSC_DIR/DATA/deepCoadd-results/HSC-I``.
2180 ConfigClass = CompareWarpAssembleCoaddConfig
2181 _DefaultName =
"compareWarpAssembleCoadd"
2184 AssembleCoaddTask.__init__(self, *args, **kwargs)
2185 self.makeSubtask(
"assembleStaticSkyModel")
2186 detectionSchema = afwTable.SourceTable.makeMinimalSchema()
2187 self.makeSubtask(
"detect", schema=detectionSchema)
2188 if self.config.doPreserveContainedBySource:
2189 self.makeSubtask(
"detectTemplate", schema=afwTable.SourceTable.makeMinimalSchema())
2190 if self.config.doScaleWarpVariance:
2191 self.makeSubtask(
"scaleWarpVariance")
2192 if self.config.doFilterMorphological:
2193 self.makeSubtask(
"maskStreaks")
2195 @utils.inheritDoc(AssembleCoaddTask)
2198 Generate a templateCoadd to use as a naive model of static sky to
2199 subtract from PSF-Matched warps.
2203 result : `lsst.pipe.base.Struct`
2204 Result struct with components:
2206 - ``templateCoadd`` : coadded exposure (``lsst.afw.image.Exposure``)
2207 - ``nImage`` : N Image (``lsst.afw.image.Image``)
2210 staticSkyModelInputRefs = copy.deepcopy(inputRefs)
2211 staticSkyModelInputRefs.inputWarps = inputRefs.psfMatchedWarps
2215 staticSkyModelOutputRefs = copy.deepcopy(outputRefs)
2216 if self.config.assembleStaticSkyModel.doWrite:
2217 staticSkyModelOutputRefs.coaddExposure = staticSkyModelOutputRefs.templateCoadd
2220 del outputRefs.templateCoadd
2221 del staticSkyModelOutputRefs.templateCoadd
2224 if 'nImage' in staticSkyModelOutputRefs.keys():
2225 del staticSkyModelOutputRefs.nImage
2227 templateCoadd = self.assembleStaticSkyModel.runQuantum(butlerQC, staticSkyModelInputRefs,
2228 staticSkyModelOutputRefs)
2229 if templateCoadd
is None:
2230 raise RuntimeError(self.
_noTemplateMessage_noTemplateMessage(self.assembleStaticSkyModel.warpType))
2232 return pipeBase.Struct(templateCoadd=templateCoadd.coaddExposure,
2233 nImage=templateCoadd.nImage,
2234 warpRefList=templateCoadd.warpRefList,
2235 imageScalerList=templateCoadd.imageScalerList,
2236 weightList=templateCoadd.weightList)
2238 @utils.inheritDoc(AssembleCoaddTask)
2241 Generate a templateCoadd to use as a naive model of static sky to
2242 subtract from PSF-Matched warps.
2246 result : `lsst.pipe.base.Struct`
2247 Result struct with components:
2249 - ``templateCoadd``: coadded exposure (``lsst.afw.image.Exposure``)
2250 - ``nImage``: N Image (``lsst.afw.image.Image``)
2252 templateCoadd = self.assembleStaticSkyModel.runDataRef(dataRef, selectDataList, warpRefList)
2253 if templateCoadd
is None:
2254 raise RuntimeError(self.
_noTemplateMessage_noTemplateMessage(self.assembleStaticSkyModel.warpType))
2256 return pipeBase.Struct(templateCoadd=templateCoadd.coaddExposure,
2257 nImage=templateCoadd.nImage,
2258 warpRefList=templateCoadd.warpRefList,
2259 imageScalerList=templateCoadd.imageScalerList,
2260 weightList=templateCoadd.weightList)
2262 def _noTemplateMessage(self, warpType):
2263 warpName = (warpType[0].upper() + warpType[1:])
2264 message =
"""No %(warpName)s warps were found to build the template coadd which is
2265 required to run CompareWarpAssembleCoaddTask. To continue assembling this type of coadd,
2266 first either rerun makeCoaddTempExp with config.make%(warpName)s=True or
2267 coaddDriver with config.makeCoadTempExp.make%(warpName)s=True, before assembleCoadd.
2269 Alternatively, to use another algorithm with existing warps, retarget the CoaddDriverConfig to
2270 another algorithm like:
2272 from lsst.pipe.tasks.assembleCoadd import SafeClipAssembleCoaddTask
2273 config.assemble.retarget(SafeClipAssembleCoaddTask)
2274 """ % {
"warpName": warpName}
2277 @utils.inheritDoc(AssembleCoaddTask)
2279 def run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
2280 supplementaryData, *args, **kwargs):
2281 """Assemble the coadd.
2283 Find artifacts and apply them to the warps' masks creating a list of
2284 alternative masks with a new "CLIPPED" plane and updated "NO_DATA"
2285 plane. Then pass these alternative masks to the base class's `run`
2288 The input parameters ``supplementaryData`` is a `lsst.pipe.base.Struct`
2289 that must contain a ``templateCoadd`` that serves as the
2290 model of the static sky.
2296 dataIds = [ref.dataId
for ref
in tempExpRefList]
2297 psfMatchedDataIds = [ref.dataId
for ref
in supplementaryData.warpRefList]
2299 if dataIds != psfMatchedDataIds:
2300 self.log.
info(
"Reordering and or/padding PSF-matched visit input list")
2301 supplementaryData.warpRefList =
reorderAndPadList(supplementaryData.warpRefList,
2302 psfMatchedDataIds, dataIds)
2303 supplementaryData.imageScalerList =
reorderAndPadList(supplementaryData.imageScalerList,
2304 psfMatchedDataIds, dataIds)
2307 spanSetMaskList = self.
findArtifactsfindArtifacts(supplementaryData.templateCoadd,
2308 supplementaryData.warpRefList,
2309 supplementaryData.imageScalerList)
2311 badMaskPlanes = self.config.badMaskPlanes[:]
2312 badMaskPlanes.append(
"CLIPPED")
2313 badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes)
2315 result = AssembleCoaddTask.run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
2316 spanSetMaskList, mask=badPixelMask)
2320 self.
applyAltEdgeMaskapplyAltEdgeMask(result.coaddExposure.maskedImage.mask, spanSetMaskList)
2324 """Propagate alt EDGE mask to SENSOR_EDGE AND INEXACT_PSF planes.
2328 mask : `lsst.afw.image.Mask`
2330 altMaskList : `list`
2331 List of Dicts containing ``spanSet`` lists.
2332 Each element contains the new mask plane name (e.g. "CLIPPED
2333 and/or "NO_DATA") as the key, and list of ``SpanSets`` to apply to
2336 maskValue = mask.getPlaneBitMask([
"SENSOR_EDGE",
"INEXACT_PSF"])
2337 for visitMask
in altMaskList:
2338 if "EDGE" in visitMask:
2339 for spanSet
in visitMask[
'EDGE']:
2340 spanSet.clippedTo(mask.getBBox()).setMask(mask, maskValue)
2345 Loop through warps twice. The first loop builds a map with the count
2346 of how many epochs each pixel deviates from the templateCoadd by more
2347 than ``config.chiThreshold`` sigma. The second loop takes each
2348 difference image and filters the artifacts detected in each using
2349 count map to filter out variable sources and sources that are
2350 difficult to subtract cleanly.
2354 templateCoadd : `lsst.afw.image.Exposure`
2355 Exposure to serve as model of static sky.
2356 tempExpRefList : `list`
2357 List of data references to warps.
2358 imageScalerList : `list`
2359 List of image scalers.
2364 List of dicts containing information about CLIPPED
2365 (i.e., artifacts), NO_DATA, and EDGE pixels.
2368 self.log.
debug(
"Generating Count Image, and mask lists.")
2369 coaddBBox = templateCoadd.getBBox()
2370 slateIm = afwImage.ImageU(coaddBBox)
2371 epochCountImage = afwImage.ImageU(coaddBBox)
2372 nImage = afwImage.ImageU(coaddBBox)
2373 spanSetArtifactList = []
2374 spanSetNoDataMaskList = []
2375 spanSetEdgeList = []
2376 spanSetBadMorphoList = []
2377 badPixelMask = self.getBadPixelMask()
2380 templateCoadd.mask.clearAllMaskPlanes()
2382 if self.config.doPreserveContainedBySource:
2383 templateFootprints = self.detectTemplate.detectFootprints(templateCoadd)
2385 templateFootprints =
None
2387 for warpRef, imageScaler
in zip(tempExpRefList, imageScalerList):
2389 if warpDiffExp
is not None:
2391 nImage.array += (numpy.isfinite(warpDiffExp.image.array)
2392 * ((warpDiffExp.mask.array & badPixelMask) == 0)).astype(numpy.uint16)
2393 fpSet = self.detect.detectFootprints(warpDiffExp, doSmooth=
False, clearMask=
True)
2394 fpSet.positive.merge(fpSet.negative)
2395 footprints = fpSet.positive
2397 spanSetList = [footprint.spans
for footprint
in footprints.getFootprints()]
2400 if self.config.doPrefilterArtifacts:
2404 self.detect.clearMask(warpDiffExp.mask)
2405 for spans
in spanSetList:
2406 spans.setImage(slateIm, 1, doClip=
True)
2407 spans.setMask(warpDiffExp.mask, warpDiffExp.mask.getPlaneBitMask(
"DETECTED"))
2408 epochCountImage += slateIm
2410 if self.config.doFilterMorphological:
2411 maskName = self.config.streakMaskName
2412 _ = self.maskStreaks.
run(warpDiffExp)
2413 streakMask = warpDiffExp.mask
2414 spanSetStreak = afwGeom.SpanSet.fromMask(streakMask,
2415 streakMask.getPlaneBitMask(maskName)).split()
2421 nans = numpy.where(numpy.isnan(warpDiffExp.maskedImage.image.array), 1, 0)
2423 nansMask.setXY0(warpDiffExp.getXY0())
2424 edgeMask = warpDiffExp.mask
2425 spanSetEdgeMask = afwGeom.SpanSet.fromMask(edgeMask,
2426 edgeMask.getPlaneBitMask(
"EDGE")).split()
2430 nansMask = afwImage.MaskX(coaddBBox, 1)
2432 spanSetEdgeMask = []
2435 spanSetNoDataMask = afwGeom.SpanSet.fromMask(nansMask).split()
2437 spanSetNoDataMaskList.append(spanSetNoDataMask)
2438 spanSetArtifactList.append(spanSetList)
2439 spanSetEdgeList.append(spanSetEdgeMask)
2440 if self.config.doFilterMorphological:
2441 spanSetBadMorphoList.append(spanSetStreak)
2444 path = self.
_dataRef2DebugPath_dataRef2DebugPath(
"epochCountIm", tempExpRefList[0], coaddLevel=
True)
2445 epochCountImage.writeFits(path)
2447 for i, spanSetList
in enumerate(spanSetArtifactList):
2449 filteredSpanSetList = self.
filterArtifactsfilterArtifacts(spanSetList, epochCountImage, nImage,
2451 spanSetArtifactList[i] = filteredSpanSetList
2452 if self.config.doFilterMorphological:
2453 spanSetArtifactList[i] += spanSetBadMorphoList[i]
2456 for artifacts, noData, edge
in zip(spanSetArtifactList, spanSetNoDataMaskList, spanSetEdgeList):
2457 altMasks.append({
'CLIPPED': artifacts,
2463 """Remove artifact candidates covered by bad mask plane.
2465 Any future editing of the candidate list that does not depend on
2466 temporal information should go in this method.
2470 spanSetList : `list`
2471 List of SpanSets representing artifact candidates.
2472 exp : `lsst.afw.image.Exposure`
2473 Exposure containing mask planes used to prefilter.
2477 returnSpanSetList : `list`
2478 List of SpanSets with artifacts.
2480 badPixelMask = exp.mask.getPlaneBitMask(self.config.prefilterArtifactsMaskPlanes)
2481 goodArr = (exp.mask.array & badPixelMask) == 0
2482 returnSpanSetList = []
2483 bbox = exp.getBBox()
2484 x0, y0 = exp.getXY0()
2485 for i, span
in enumerate(spanSetList):
2486 y, x = span.clippedTo(bbox).indices()
2487 yIndexLocal = numpy.array(y) - y0
2488 xIndexLocal = numpy.array(x) - x0
2489 goodRatio = numpy.count_nonzero(goodArr[yIndexLocal, xIndexLocal])/span.getArea()
2490 if goodRatio > self.config.prefilterArtifactsRatio:
2491 returnSpanSetList.append(span)
2492 return returnSpanSetList
2494 def filterArtifacts(self, spanSetList, epochCountImage, nImage, footprintsToExclude=None):
2495 """Filter artifact candidates.
2499 spanSetList : `list`
2500 List of SpanSets representing artifact candidates.
2501 epochCountImage : `lsst.afw.image.Image`
2502 Image of accumulated number of warpDiff detections.
2503 nImage : `lsst.afw.image.Image`
2504 Image of the accumulated number of total epochs contributing.
2508 maskSpanSetList : `list`
2509 List of SpanSets with artifacts.
2512 maskSpanSetList = []
2513 x0, y0 = epochCountImage.getXY0()
2514 for i, span
in enumerate(spanSetList):
2515 y, x = span.indices()
2516 yIdxLocal = [y1 - y0
for y1
in y]
2517 xIdxLocal = [x1 - x0
for x1
in x]
2518 outlierN = epochCountImage.array[yIdxLocal, xIdxLocal]
2519 totalN = nImage.array[yIdxLocal, xIdxLocal]
2522 effMaxNumEpochsHighN = (self.config.maxNumEpochs
2523 + self.config.maxFractionEpochsHigh*numpy.mean(totalN))
2524 effMaxNumEpochsLowN = self.config.maxFractionEpochsLow * numpy.mean(totalN)
2525 effectiveMaxNumEpochs = int(
min(effMaxNumEpochsLowN, effMaxNumEpochsHighN))
2526 nPixelsBelowThreshold = numpy.count_nonzero((outlierN > 0)
2527 & (outlierN <= effectiveMaxNumEpochs))
2528 percentBelowThreshold = nPixelsBelowThreshold / len(outlierN)
2529 if percentBelowThreshold > self.config.spatialThreshold:
2530 maskSpanSetList.append(span)
2532 if self.config.doPreserveContainedBySource
and footprintsToExclude
is not None:
2534 filteredMaskSpanSetList = []
2535 for span
in maskSpanSetList:
2537 for footprint
in footprintsToExclude.positive.getFootprints():
2538 if footprint.spans.contains(span):
2542 filteredMaskSpanSetList.append(span)
2543 maskSpanSetList = filteredMaskSpanSetList
2545 return maskSpanSetList
2547 def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd):
2548 """Fetch a warp from the butler and return a warpDiff.
2552 warpRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
2553 Butler dataRef for the warp.
2554 imageScaler : `lsst.pipe.tasks.scaleZeroPoint.ImageScaler`
2555 An image scaler object.
2556 templateCoadd : `lsst.afw.image.Exposure`
2557 Exposure to be substracted from the scaled warp.
2561 warp : `lsst.afw.image.Exposure`
2562 Exposure of the image difference between the warp and template.
2570 warpName = self.getTempExpDatasetName(
'psfMatched')
2571 if not isinstance(warpRef, DeferredDatasetHandle):
2572 if not warpRef.datasetExists(warpName):
2573 self.log.
warning(
"Could not find %s %s; skipping it", warpName, warpRef.dataId)
2575 warp = warpRef.get(datasetType=warpName, immediate=
True)
2577 imageScaler.scaleMaskedImage(warp.getMaskedImage())
2578 mi = warp.getMaskedImage()
2579 if self.config.doScaleWarpVariance:
2581 self.scaleWarpVariance.
run(mi)
2582 except Exception
as exc:
2583 self.log.
warning(
"Unable to rescale variance of warp (%s); leaving it as-is", exc)
2584 mi -= templateCoadd.getMaskedImage()
2587 def _dataRef2DebugPath(self, prefix, warpRef, coaddLevel=False):
2588 """Return a path to which to write debugging output.
2590 Creates a hyphen-delimited string of dataId values for simple filenames.
2595 Prefix for filename.
2596 warpRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
2597 Butler dataRef to make the path from.
2598 coaddLevel : `bool`, optional.
2599 If True, include only coadd-level keys (e.g., 'tract', 'patch',
2600 'filter', but no 'visit').
2605 Path for debugging output.
2608 keys = warpRef.getButler().getKeys(self.getCoaddDatasetName(self.warpType))
2610 keys = warpRef.dataId.keys()
2611 keyList = sorted(keys, reverse=
True)
2613 filename =
"%s-%s.fits" % (prefix,
'-'.join([str(warpRef.dataId[k])
for k
in keyList]))
2614 return os.path.join(directory, filename)
std::vector< SchemaItem< Flag > > * items
A Threshold is used to pass a threshold value to detection algorithms.
A compact representation of a collection of pixels.
A group of labels for a filter in an exposure or coadd.
Represent a 2-dimensional array of bitmask pixels.
Pass parameters to a Statistics object.
A floating-point coordinate rectangle geometry.
An integer coordinate rectangle.
Reports invalid arguments.
def makeDataRefList(self, namespace)
def __init__(self, *config=None)
def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs)
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, supplementaryData, *args, **kwargs)
def prefilterArtifacts(self, spanSetList, exp)
def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList)
def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd)
def _dataRef2DebugPath(self, prefix, warpRef, coaddLevel=False)
def applyAltEdgeMask(self, mask, altMaskList)
def makeSupplementaryData(self, dataRef, selectDataList=None, warpRefList=None)
def _noTemplateMessage(self, warpType)
def filterArtifacts(self, spanSetList, epochCountImage, nImage, footprintsToExclude=None)
def __init__(self, *args, **kwargs)
def buildDifferenceImage(self, skyInfo, tempExpRefList, imageScalerList, weightList)
def detectClipBig(self, clipList, clipFootprints, clipIndices, detectionFootprints, maskClipValue, maskDetValue, coaddBBox)
def __init__(self, *args, **kwargs)
def detectClip(self, exp, tempExpRefList)
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, *args, **kwargs)
Base class for coaddition.
daf::base::PropertyList * list
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
std::shared_ptr< lsst::afw::image::Image< PixelT > > statisticsStack(std::vector< std::shared_ptr< lsst::afw::image::Image< PixelT >>> &images, Property flags, StatisticsControl const &sctrl=StatisticsControl(), std::vector< lsst::afw::image::VariancePixel > const &wvector=std::vector< lsst::afw::image::VariancePixel >(0))
A function to compute some statistics of a stack of Images.
void setCoaddEdgeBits(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > &coaddMask, lsst::afw::image::Image< WeightPixelT > const &weightMap)
set edge bits of coadd mask based on weight map
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
def reorderAndPadList(inputList, inputKeys, outputKeys, padWith=None)
def prepareStats(self, mask=None)
def readBrightObjectMasks(self, dataRef)
def makeSupplementaryData(self, dataRef, selectDataList=None, warpRefList=None)
def countMaskFromFootprint(mask, footprint, bitmask, ignoreMask)
def assembleOnlineMeanCoadd(self, coaddExposure, tempExpRefList, imageScalerList, weightList, altMaskList, statsCtrl, nImage=None)
def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs)
def applyAltMaskPlanes(self, mask, altMaskSpans)
def shrinkValidPolygons(self, coaddInputs)
def setBrightObjectMasks(self, exposure, brightObjectMasks, dataId=None)
def getTempExpRefList(self, patchRef, calExpRefList)
def assembleMetadata(self, coaddExposure, tempExpRefList, weightList)
def removeMaskPlanes(self, maskedImage)
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)
def filterWarps(self, inputs, goodVisits)
def prepareInputs(self, refList)
def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList, weightList, altMaskList, statsFlags, statsCtrl, nImage=None)
def setInexactPsf(self, mask)
def processResults(self, coaddExposure, brightObjectMasks=None, dataId=None)
def makeCoaddSuffix(warpType="direct")
def makeSkyInfo(skyMap, tractId, patchId)
def getGroupDataRef(butler, datasetType, groupTuple, keys)
def groupPatchExposures(patchDataRef, calexpDataRefList, coaddDatasetType="deepCoadd", tempExpDatasetType="deepCoadd_directWarp")