25 from scipy
import ndimage
30 from lsst.daf.butler
import DeferredDatasetHandle
37 from lsst.utils.timer
import timeMethod
38 from .assembleCoadd
import (AssembleCoaddConnections,
40 CompareWarpAssembleCoaddConfig,
41 CompareWarpAssembleCoaddTask)
42 from .coaddBase
import makeSkyInfo
43 from .measurePsf
import MeasurePsfTask
45 __all__ = [
"DcrAssembleCoaddConnections",
"DcrAssembleCoaddTask",
"DcrAssembleCoaddConfig"]
49 dimensions=(
"tract",
"patch",
"band",
"skymap"),
50 defaultTemplates={
"inputWarpName":
"deep",
51 "inputCoaddName":
"deep",
52 "outputCoaddName":
"dcr",
56 inputWarps = pipeBase.connectionTypes.Input(
57 doc=(
"Input list of warps to be assembled i.e. stacked."
58 "Note that this will often be different than the inputCoaddName."
59 "WarpType (e.g. direct, psfMatched) is controlled by the warpType config parameter"),
60 name=
"{inputWarpName}Coadd_{warpType}Warp",
61 storageClass=
"ExposureF",
62 dimensions=(
"tract",
"patch",
"skymap",
"visit",
"instrument"),
66 templateExposure = pipeBase.connectionTypes.Input(
67 doc=
"Input coadded exposure, produced by previous call to AssembleCoadd",
68 name=
"{fakesType}{inputCoaddName}Coadd{warpTypeSuffix}",
69 storageClass=
"ExposureF",
70 dimensions=(
"tract",
"patch",
"skymap",
"band"),
72 dcrCoadds = pipeBase.connectionTypes.Output(
73 doc=
"Output coadded exposure, produced by stacking input warps",
74 name=
"{fakesType}{outputCoaddName}Coadd{warpTypeSuffix}",
75 storageClass=
"ExposureF",
76 dimensions=(
"tract",
"patch",
"skymap",
"band",
"subfilter"),
79 dcrNImages = pipeBase.connectionTypes.Output(
80 doc=
"Output image of number of input images per pixel",
81 name=
"{outputCoaddName}Coadd_nImage",
82 storageClass=
"ImageU",
83 dimensions=(
"tract",
"patch",
"skymap",
"band",
"subfilter"),
87 def __init__(self, *, config=None):
88 super().__init__(config=config)
89 if not config.doWrite:
90 self.outputs.remove(
"dcrCoadds")
91 if not config.doNImage:
92 self.outputs.remove(
"dcrNImages")
94 self.outputs.remove(
"coaddExposure")
95 self.outputs.remove(
"nImage")
99 pipelineConnections=DcrAssembleCoaddConnections):
100 dcrNumSubfilters = pexConfig.Field(
102 doc=
"Number of sub-filters to forward model chromatic effects to fit the supplied exposures.",
105 maxNumIter = pexConfig.Field(
108 doc=
"Maximum number of iterations of forward modeling.",
111 minNumIter = pexConfig.Field(
114 doc=
"Minimum number of iterations of forward modeling.",
117 convergenceThreshold = pexConfig.Field(
119 doc=
"Target relative change in convergence between iterations of forward modeling.",
122 useConvergence = pexConfig.Field(
124 doc=
"Use convergence test as a forward modeling end condition?"
125 "If not set, skips calculating convergence and runs for ``maxNumIter`` iterations",
128 baseGain = pexConfig.Field(
131 doc=
"Relative weight to give the new solution vs. the last solution when updating the model."
132 "A value of 1.0 gives equal weight to both solutions."
133 "Small values imply slower convergence of the solution, but can "
134 "help prevent overshooting and failures in the fit."
135 "If ``baseGain`` is None, a conservative gain "
136 "will be calculated from the number of subfilters. ",
139 useProgressiveGain = pexConfig.Field(
141 doc=
"Use a gain that slowly increases above ``baseGain`` to accelerate convergence? "
142 "When calculating the next gain, we use up to 5 previous gains and convergence values."
143 "Can be set to False to force the model to change at the rate of ``baseGain``. ",
146 doAirmassWeight = pexConfig.Field(
148 doc=
"Weight exposures by airmass? Useful if there are relatively few high-airmass observations.",
151 modelWeightsWidth = pexConfig.Field(
153 doc=
"Width of the region around detected sources to include in the DcrModel.",
156 useModelWeights = pexConfig.Field(
158 doc=
"Width of the region around detected sources to include in the DcrModel.",
161 splitSubfilters = pexConfig.Field(
163 doc=
"Calculate DCR for two evenly-spaced wavelengths in each subfilter."
164 "Instead of at the midpoint",
167 splitThreshold = pexConfig.Field(
169 doc=
"Minimum DCR difference within a subfilter to use ``splitSubfilters``, in pixels."
170 "Set to 0 to always split the subfilters.",
173 regularizeModelIterations = pexConfig.Field(
175 doc=
"Maximum relative change of the model allowed between iterations."
176 "Set to zero to disable.",
179 regularizeModelFrequency = pexConfig.Field(
181 doc=
"Maximum relative change of the model allowed between subfilters."
182 "Set to zero to disable.",
185 convergenceMaskPlanes = pexConfig.ListField(
187 default=[
"DETECTED"],
188 doc=
"Mask planes to use to calculate convergence."
190 regularizationWidth = pexConfig.Field(
193 doc=
"Minimum radius of a region to include in regularization, in pixels."
195 imageInterpOrder = pexConfig.Field(
197 doc=
"The order of the spline interpolation used to shift the image plane.",
200 accelerateModel = pexConfig.Field(
202 doc=
"Factor to amplify the differences between model planes by to speed convergence.",
205 doCalculatePsf = pexConfig.Field(
207 doc=
"Set to detect stars and recalculate the PSF from the final coadd."
208 "Otherwise the PSF is estimated from a selection of the best input exposures",
211 detectPsfSources = pexConfig.ConfigurableField(
212 target=measAlg.SourceDetectionTask,
213 doc=
"Task to detect sources for PSF measurement, if ``doCalculatePsf`` is set.",
215 measurePsfSources = pexConfig.ConfigurableField(
216 target=SingleFrameMeasurementTask,
217 doc=
"Task to measure sources for PSF measurement, if ``doCalculatePsf`` is set."
219 measurePsf = pexConfig.ConfigurableField(
220 target=MeasurePsfTask,
221 doc=
"Task to measure the PSF of the coadd, if ``doCalculatePsf`` is set.",
223 effectiveWavelength = pexConfig.Field(
224 doc=
"Effective wavelength of the filter, in nm."
225 "Required if transmission curves aren't used."
226 "Support for using transmission curves is to be added in DM-13668.",
229 bandwidth = pexConfig.Field(
230 doc=
"Bandwidth of the physical filter, in nm."
231 "Required if transmission curves aren't used."
232 "Support for using transmission curves is to be added in DM-13668.",
237 CompareWarpAssembleCoaddConfig.setDefaults(self)
238 self.assembleStaticSkyModel.retarget(CompareWarpAssembleCoaddTask)
240 self.assembleStaticSkyModel.warpType = self.warpType
242 self.assembleStaticSkyModel.doNImage =
False
243 self.assembleStaticSkyModel.doWrite =
False
244 self.detectPsfSources.returnOriginalFootprints =
False
245 self.detectPsfSources.thresholdPolarity =
"positive"
247 self.detectPsfSources.thresholdValue = 50
248 self.detectPsfSources.nSigmaToGrow = 2
250 self.detectPsfSources.minPixels = 25
252 self.detectPsfSources.thresholdType =
"pixel_stdev"
255 self.measurePsf.starSelector[
"objectSize"].doFluxLimit =
False
259 """Assemble DCR coadded images from a set of warps.
264 The number of pixels to grow each subregion by to allow for DCR.
268 As with AssembleCoaddTask, we want to assemble a coadded image from a set of
269 Warps (also called coadded temporary exposures), including the effects of
270 Differential Chromatic Refraction (DCR).
271 For full details of the mathematics and algorithm, please see
272 DMTN-037: DCR-matched template generation (https://dmtn-037.lsst.io).
274 This Task produces a DCR-corrected deepCoadd, as well as a dcrCoadd for
275 each subfilter used in the iterative calculation.
276 It begins by dividing the bandpass-defining filter into N equal bandwidth
277 "subfilters", and divides the flux in each pixel from an initial coadd
278 equally into each as a "dcrModel". Because the airmass and parallactic
279 angle of each individual exposure is known, we can calculate the shift
280 relative to the center of the band in each subfilter due to DCR. For each
281 exposure we apply this shift as a linear transformation to the dcrModels
282 and stack the results to produce a DCR-matched exposure. The matched
283 exposures are subtracted from the input exposures to produce a set of
284 residual images, and these residuals are reverse shifted for each
285 exposures' subfilters and stacked. The shifted and stacked residuals are
286 added to the dcrModels to produce a new estimate of the flux in each pixel
287 within each subfilter. The dcrModels are solved for iteratively, which
288 continues until the solution from a new iteration improves by less than
289 a set percentage, or a maximum number of iterations is reached.
290 Two forms of regularization are employed to reduce unphysical results.
291 First, the new solution is averaged with the solution from the previous
292 iteration, which mitigates oscillating solutions where the model
293 overshoots with alternating very high and low values.
294 Second, a common degeneracy when the data have a limited range of airmass or
295 parallactic angle values is for one subfilter to be fit with very low or
296 negative values, while another subfilter is fit with very high values. This
297 typically appears in the form of holes next to sources in one subfilter,
298 and corresponding extended wings in another. Because each subfilter has
299 a narrow bandwidth we assume that physical sources that are above the noise
300 level will not vary in flux by more than a factor of `frequencyClampFactor`
301 between subfilters, and pixels that have flux deviations larger than that
302 factor will have the excess flux distributed evenly among all subfilters.
303 If `splitSubfilters` is set, then each subfilter will be further sub-
304 divided during the forward modeling step (only). This approximates using
305 a higher number of subfilters that may be necessary for high airmass
306 observations, but does not increase the number of free parameters in the
307 fit. This is needed when there are high airmass observations which would
308 otherwise have significant DCR even within a subfilter. Because calculating
309 the shifted images takes most of the time, splitting the subfilters is
310 turned off by way of the `splitThreshold` option for low-airmass
311 observations that do not suffer from DCR within a subfilter.
314 ConfigClass = DcrAssembleCoaddConfig
315 _DefaultName =
"dcrAssembleCoadd"
317 def __init__(self, *args, **kwargs):
318 super().__init__(*args, **kwargs)
319 if self.config.doCalculatePsf:
320 self.schema = afwTable.SourceTable.makeMinimalSchema()
321 self.makeSubtask(
"detectPsfSources", schema=self.schema)
322 self.makeSubtask(
"measurePsfSources", schema=self.schema)
323 self.makeSubtask(
"measurePsf", schema=self.schema)
325 @utils.inheritDoc(pipeBase.PipelineTask)
326 def runQuantum(self, butlerQC, inputRefs, outputRefs):
331 Assemble a coadd from a set of Warps.
333 PipelineTask (Gen3) entry point to Coadd a set of Warps.
334 Analogous to `runDataRef`, it prepares all the data products to be
335 passed to `run`, and processes the results before returning a struct
336 of results to be written out. AssembleCoadd cannot fit all Warps in memory.
337 Therefore, its inputs are accessed subregion by subregion
338 by the Gen3 `DeferredDatasetHandle` that is analagous to the Gen2
339 `lsst.daf.persistence.ButlerDataRef`. Any updates to this method should
340 correspond to an update in `runDataRef` while both entry points
343 inputData = butlerQC.get(inputRefs)
347 skyMap = inputData[
"skyMap"]
348 outputDataId = butlerQC.quantum.dataId
351 tractId=outputDataId[
'tract'],
352 patchId=outputDataId[
'patch'])
356 warpRefList = inputData[
'inputWarps']
358 inputs = self.prepareInputs(warpRefList)
359 self.log.
info(
"Found %d %s", len(inputs.tempExpRefList),
360 self.getTempExpDatasetName(self.warpType))
361 if len(inputs.tempExpRefList) == 0:
362 self.log.
warning(
"No coadd temporary exposures found")
365 supplementaryData = self.makeSupplementaryDataGen3(butlerQC, inputRefs, outputRefs)
366 retStruct = self.run(inputData[
'skyInfo'], inputs.tempExpRefList, inputs.imageScalerList,
367 inputs.weightList, supplementaryData=supplementaryData)
369 inputData.setdefault(
'brightObjectMask',
None)
370 for subfilter
in range(self.config.dcrNumSubfilters):
372 retStruct.dcrCoadds[subfilter].setPsf(retStruct.coaddExposure.getPsf())
373 self.processResults(retStruct.dcrCoadds[subfilter], inputData[
'brightObjectMask'], outputDataId)
375 if self.config.doWrite:
376 butlerQC.put(retStruct, outputRefs)
380 def runDataRef(self, dataRef, selectDataList=None, warpRefList=None):
381 """Assemble a coadd from a set of warps.
383 Coadd a set of Warps. Compute weights to be applied to each Warp and
384 find scalings to match the photometric zeropoint to a reference Warp.
385 Assemble the Warps using run method.
386 Forward model chromatic effects across multiple subfilters,
387 and subtract from the input Warps to build sets of residuals.
388 Use the residuals to construct a new ``DcrModel`` for each subfilter,
389 and iterate until the model converges.
390 Interpolate over NaNs and optionally write the coadd to disk.
391 Return the coadded exposure.
395 dataRef : `lsst.daf.persistence.ButlerDataRef`
396 Data reference defining the patch for coaddition and the
398 selectDataList : `list` of `lsst.daf.persistence.ButlerDataRef`
399 List of data references to warps. Data to be coadded will be
400 selected from this list based on overlap with the patch defined by
405 results : `lsst.pipe.base.Struct`
406 The Struct contains the following fields:
408 - ``coaddExposure``: coadded exposure (`lsst.afw.image.Exposure`)
409 - ``nImage``: exposure count image (`lsst.afw.image.ImageU`)
410 - ``dcrCoadds``: `list` of coadded exposures for each subfilter
411 - ``dcrNImages``: `list` of exposure count images for each subfilter
413 if (selectDataList
is None and warpRefList
is None)
or (selectDataList
and warpRefList):
414 raise RuntimeError(
"runDataRef must be supplied either a selectDataList or warpRefList")
416 skyInfo = self.getSkyInfo(dataRef)
417 if warpRefList
is None:
418 calExpRefList = self.selectExposures(dataRef, skyInfo, selectDataList=selectDataList)
419 if len(calExpRefList) == 0:
420 self.log.
warning(
"No exposures to coadd")
422 self.log.
info(
"Coadding %d exposures", len(calExpRefList))
424 warpRefList = self.getTempExpRefList(dataRef, calExpRefList)
426 inputData = self.prepareInputs(warpRefList)
427 self.log.
info(
"Found %d %s", len(inputData.tempExpRefList),
428 self.getTempExpDatasetName(self.warpType))
429 if len(inputData.tempExpRefList) == 0:
430 self.log.
warning(
"No coadd temporary exposures found")
433 supplementaryData = self.makeSupplementaryData(dataRef, warpRefList=inputData.tempExpRefList)
435 results = self.run(skyInfo, inputData.tempExpRefList, inputData.imageScalerList,
436 inputData.weightList, supplementaryData=supplementaryData)
438 self.log.
warning(
"Could not construct DcrModel for patch %s: no data to coadd.",
439 skyInfo.patchInfo.getIndex())
442 if self.config.doCalculatePsf:
443 self.measureCoaddPsf(results.coaddExposure)
444 brightObjects = self.readBrightObjectMasks(dataRef)
if self.config.doMaskBrightObjects
else None
445 for subfilter
in range(self.config.dcrNumSubfilters):
447 results.dcrCoadds[subfilter].setPsf(results.coaddExposure.getPsf())
448 self.processResults(results.dcrCoadds[subfilter],
449 brightObjectMasks=brightObjects, dataId=dataRef.dataId)
450 if self.config.doWrite:
451 self.log.
info(
"Persisting dcrCoadd")
452 dataRef.put(results.dcrCoadds[subfilter],
"dcrCoadd", subfilter=subfilter,
453 numSubfilters=self.config.dcrNumSubfilters)
454 if self.config.doNImage
and results.dcrNImages
is not None:
455 dataRef.put(results.dcrNImages[subfilter],
"dcrCoadd_nImage", subfilter=subfilter,
456 numSubfilters=self.config.dcrNumSubfilters)
460 @utils.inheritDoc(AssembleCoaddTask)
462 """Load the previously-generated template coadd.
464 This can be removed entirely once we no longer support the Gen 2 butler.
468 templateCoadd : `lsst.pipe.base.Struct`
469 Result struct with components:
471 - ``templateCoadd``: coadded exposure (`lsst.afw.image.ExposureF`)
473 templateCoadd = butlerQC.get(inputRefs.templateExposure)
475 return pipeBase.Struct(templateCoadd=templateCoadd)
477 def measureCoaddPsf(self, coaddExposure):
478 """Detect sources on the coadd exposure and measure the final PSF.
482 coaddExposure : `lsst.afw.image.Exposure`
483 The final coadded exposure.
485 table = afwTable.SourceTable.make(self.schema)
486 detResults = self.detectPsfSources.
run(table, coaddExposure, clearMask=
False)
487 coaddSources = detResults.sources
488 self.measurePsfSources.
run(
489 measCat=coaddSources,
490 exposure=coaddExposure
497 psfResults = self.measurePsf.
run(coaddExposure, coaddSources)
498 except Exception
as e:
499 self.log.
warning(
"Unable to calculate PSF, using default coadd PSF: %s", e)
501 coaddExposure.setPsf(psfResults.psf)
503 def prepareDcrInputs(self, templateCoadd, warpRefList, weightList):
504 """Prepare the DCR coadd by iterating through the visitInfo of the input warps.
506 Sets the property ``bufferSize``.
510 templateCoadd : `lsst.afw.image.ExposureF`
511 The initial coadd exposure before accounting for DCR.
512 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
513 `lsst.daf.persistence.ButlerDataRef`
514 The data references to the input warped exposures.
515 weightList : `list` of `float`
516 The weight to give each input exposure in the coadd
517 Will be modified in place if ``doAirmassWeight`` is set.
521 dcrModels : `lsst.pipe.tasks.DcrModel`
522 Best fit model of the true sky after correcting chromatic effects.
527 If ``lambdaMin`` is missing from the Mapper class of the obs package being used.
529 sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
530 filterLabel = templateCoadd.getFilterLabel()
531 tempExpName = self.getTempExpDatasetName(self.warpType)
536 for visitNum, warpExpRef
in enumerate(warpRefList):
537 if isinstance(warpExpRef, DeferredDatasetHandle):
539 visitInfo = warpExpRef.get(component=
"visitInfo")
540 psf = warpExpRef.get(component=
"psf")
543 visitInfo = warpExpRef.get(tempExpName +
"_visitInfo")
544 psf = warpExpRef.get(tempExpName).getPsf()
545 visit = warpExpRef.dataId[
"visit"]
546 psfSize = psf.computeShape().getDeterminantRadius()*sigma2fwhm
547 airmass = visitInfo.getBoresightAirmass()
548 parallacticAngle = visitInfo.getBoresightParAngle().asDegrees()
549 airmassDict[visit] = airmass
550 angleDict[visit] = parallacticAngle
551 psfSizeDict[visit] = psfSize
552 if self.config.doAirmassWeight:
553 weightList[visitNum] *= airmass
554 dcrShifts.append(np.max(np.abs(
calculateDcr(visitInfo, templateCoadd.getWcs(),
555 self.config.effectiveWavelength,
556 self.config.bandwidth,
557 self.config.dcrNumSubfilters))))
558 self.log.
info(
"Selected airmasses:\n%s", airmassDict)
559 self.log.
info(
"Selected parallactic angles:\n%s", angleDict)
560 self.log.
info(
"Selected PSF sizes:\n%s", psfSizeDict)
561 self.bufferSize = int(np.ceil(np.max(dcrShifts)) + 1)
563 psf = self.selectCoaddPsf(templateCoadd, warpRefList)
564 except Exception
as e:
565 self.log.
warning(
"Unable to calculate restricted PSF, using default coadd PSF: %s", e)
567 psf = templateCoadd.getPsf()
568 dcrModels = DcrModel.fromImage(templateCoadd.maskedImage,
569 self.config.dcrNumSubfilters,
570 effectiveWavelength=self.config.effectiveWavelength,
571 bandwidth=self.config.bandwidth,
572 filterLabel=filterLabel,
577 def run(self, skyInfo, warpRefList, imageScalerList, weightList,
578 supplementaryData=None):
579 """Assemble the coadd.
581 Requires additional inputs Struct ``supplementaryData`` to contain a
582 ``templateCoadd`` that serves as the model of the static sky.
584 Find artifacts and apply them to the warps' masks creating a list of
585 alternative masks with a new "CLIPPED" plane and updated "NO_DATA" plane
586 Then pass these alternative masks to the base class's assemble method.
588 Divide the ``templateCoadd`` evenly between each subfilter of a
589 ``DcrModel`` as the starting best estimate of the true wavelength-
590 dependent sky. Forward model the ``DcrModel`` using the known
591 chromatic effects in each subfilter and calculate a convergence metric
592 based on how well the modeled template matches the input warps. If
593 the convergence has not yet reached the desired threshold, then shift
594 and stack the residual images to build a new ``DcrModel``. Apply
595 conditioning to prevent oscillating solutions between iterations or
598 Once the ``DcrModel`` reaches convergence or the maximum number of
599 iterations has been reached, fill the metadata for each subfilter
600 image and make them proper ``coaddExposure``s.
604 skyInfo : `lsst.pipe.base.Struct`
605 Patch geometry information, from getSkyInfo
606 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
607 `lsst.daf.persistence.ButlerDataRef`
608 The data references to the input warped exposures.
609 imageScalerList : `list` of `lsst.pipe.task.ImageScaler`
610 The image scalars correct for the zero point of the exposures.
611 weightList : `list` of `float`
612 The weight to give each input exposure in the coadd
613 supplementaryData : `lsst.pipe.base.Struct`
614 Result struct returned by ``makeSupplementaryData`` with components:
616 - ``templateCoadd``: coadded exposure (`lsst.afw.image.Exposure`)
620 result : `lsst.pipe.base.Struct`
621 Result struct with components:
623 - ``coaddExposure``: coadded exposure (`lsst.afw.image.Exposure`)
624 - ``nImage``: exposure count image (`lsst.afw.image.ImageU`)
625 - ``dcrCoadds``: `list` of coadded exposures for each subfilter
626 - ``dcrNImages``: `list` of exposure count images for each subfilter
628 minNumIter = self.config.minNumIter
or self.config.dcrNumSubfilters
629 maxNumIter = self.config.maxNumIter
or self.config.dcrNumSubfilters*3
630 templateCoadd = supplementaryData.templateCoadd
631 baseMask = templateCoadd.mask.clone()
634 baseVariance = templateCoadd.variance.clone()
635 baseVariance /= self.config.dcrNumSubfilters
636 spanSetMaskList = self.findArtifacts(templateCoadd, warpRefList, imageScalerList)
638 templateCoadd.setMask(baseMask)
639 badMaskPlanes = self.config.badMaskPlanes[:]
644 badPixelMask = templateCoadd.mask.getPlaneBitMask(badMaskPlanes)
646 stats = self.prepareStats(mask=badPixelMask)
647 dcrModels = self.prepareDcrInputs(templateCoadd, warpRefList, weightList)
648 if self.config.doNImage:
649 dcrNImages, dcrWeights = self.calculateNImage(dcrModels, skyInfo.bbox, warpRefList,
650 spanSetMaskList, stats.ctrl)
651 nImage = afwImage.ImageU(skyInfo.bbox)
655 for dcrNImage
in dcrNImages:
661 nSubregions = (
ceil(skyInfo.bbox.getHeight()/subregionSize[1])
662 *
ceil(skyInfo.bbox.getWidth()/subregionSize[0]))
664 for subBBox
in self._subBBoxIter(skyInfo.bbox, subregionSize):
667 self.log.
info(
"Computing coadd over patch %s subregion %s of %s: %s",
668 skyInfo.patchInfo.getIndex(), subIter, nSubregions, subBBox)
670 dcrBBox.grow(self.bufferSize)
671 dcrBBox.clip(dcrModels.bbox)
672 modelWeights = self.calculateModelWeights(dcrModels, dcrBBox)
673 subExposures = self.loadSubExposures(dcrBBox, stats.ctrl, warpRefList,
674 imageScalerList, spanSetMaskList)
675 convergenceMetric = self.calculateConvergence(dcrModels, subExposures, subBBox,
676 warpRefList, weightList, stats.ctrl)
677 self.log.
info(
"Initial convergence : %s", convergenceMetric)
678 convergenceList = [convergenceMetric]
680 convergenceCheck = 1.
681 refImage = templateCoadd.image
682 while (convergenceCheck > self.config.convergenceThreshold
or modelIter <= minNumIter):
683 gain = self.calculateGain(convergenceList, gainList)
684 self.dcrAssembleSubregion(dcrModels, subExposures, subBBox, dcrBBox, warpRefList,
685 stats.ctrl, convergenceMetric, gain,
686 modelWeights, refImage, dcrWeights)
687 if self.config.useConvergence:
688 convergenceMetric = self.calculateConvergence(dcrModels, subExposures, subBBox,
689 warpRefList, weightList, stats.ctrl)
690 if convergenceMetric == 0:
691 self.log.
warning(
"Coadd patch %s subregion %s had convergence metric of 0.0 which is "
692 "most likely due to there being no valid data in the region.",
693 skyInfo.patchInfo.getIndex(), subIter)
695 convergenceCheck = (convergenceList[-1] - convergenceMetric)/convergenceMetric
696 if (convergenceCheck < 0) & (modelIter > minNumIter):
697 self.log.
warning(
"Coadd patch %s subregion %s diverged before reaching maximum "
698 "iterations or desired convergence improvement of %s."
700 skyInfo.patchInfo.getIndex(), subIter,
701 self.config.convergenceThreshold, convergenceCheck)
703 convergenceList.append(convergenceMetric)
704 if modelIter > maxNumIter:
705 if self.config.useConvergence:
706 self.log.
warning(
"Coadd patch %s subregion %s reached maximum iterations "
707 "before reaching desired convergence improvement of %s."
708 " Final convergence improvement: %s",
709 skyInfo.patchInfo.getIndex(), subIter,
710 self.config.convergenceThreshold, convergenceCheck)
713 if self.config.useConvergence:
714 self.log.
info(
"Iteration %s with convergence metric %s, %.4f%% improvement (gain: %.2f)",
715 modelIter, convergenceMetric, 100.*convergenceCheck, gain)
718 if self.config.useConvergence:
719 self.log.
info(
"Coadd patch %s subregion %s finished with "
720 "convergence metric %s after %s iterations",
721 skyInfo.patchInfo.getIndex(), subIter, convergenceMetric, modelIter)
723 self.log.
info(
"Coadd patch %s subregion %s finished after %s iterations",
724 skyInfo.patchInfo.getIndex(), subIter, modelIter)
725 if self.config.useConvergence
and convergenceMetric > 0:
726 self.log.
info(
"Final convergence improvement was %.4f%% overall",
727 100*(convergenceList[0] - convergenceMetric)/convergenceMetric)
729 dcrCoadds = self.fillCoadd(dcrModels, skyInfo, warpRefList, weightList,
730 calibration=self.scaleZeroPoint.getPhotoCalib(),
731 coaddInputs=templateCoadd.getInfo().getCoaddInputs(),
733 variance=baseVariance)
734 coaddExposure = self.stackCoadd(dcrCoadds)
735 return pipeBase.Struct(coaddExposure=coaddExposure, nImage=nImage,
736 dcrCoadds=dcrCoadds, dcrNImages=dcrNImages)
738 def calculateNImage(self, dcrModels, bbox, warpRefList, spanSetMaskList, statsCtrl):
739 """Calculate the number of exposures contributing to each subfilter.
743 dcrModels : `lsst.pipe.tasks.DcrModel`
744 Best fit model of the true sky after correcting chromatic effects.
745 bbox : `lsst.geom.box.Box2I`
746 Bounding box of the patch to coadd.
747 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
748 `lsst.daf.persistence.ButlerDataRef`
749 The data references to the input warped exposures.
750 spanSetMaskList : `list` of `dict` containing spanSet lists, or None
751 Each element of the `dict` contains the new mask plane name
752 (e.g. "CLIPPED and/or "NO_DATA") as the key,
753 and the list of SpanSets to apply to the mask.
754 statsCtrl : `lsst.afw.math.StatisticsControl`
755 Statistics control object for coadd
759 dcrNImages : `list` of `lsst.afw.image.ImageU`
760 List of exposure count images for each subfilter
761 dcrWeights : `list` of `lsst.afw.image.ImageF`
762 Per-pixel weights for each subfilter.
763 Equal to 1/(number of unmasked images contributing to each pixel).
765 dcrNImages = [afwImage.ImageU(bbox)
for subfilter
in range(self.config.dcrNumSubfilters)]
766 dcrWeights = [afwImage.ImageF(bbox)
for subfilter
in range(self.config.dcrNumSubfilters)]
767 tempExpName = self.getTempExpDatasetName(self.warpType)
768 for warpExpRef, altMaskSpans
in zip(warpRefList, spanSetMaskList):
769 if isinstance(warpExpRef, DeferredDatasetHandle):
771 exposure = warpExpRef.get(parameters={
'bbox': bbox})
774 exposure = warpExpRef.get(tempExpName +
"_sub", bbox=bbox)
775 visitInfo = exposure.getInfo().getVisitInfo()
776 wcs = exposure.getInfo().getWcs()
778 if altMaskSpans
is not None:
779 self.applyAltMaskPlanes(mask, altMaskSpans)
780 weightImage = np.zeros_like(exposure.image.array)
781 weightImage[(mask.array & statsCtrl.getAndMask()) == 0] = 1.
784 weightsGenerator = self.dcrResiduals(weightImage, visitInfo, wcs,
785 dcrModels.effectiveWavelength, dcrModels.bandwidth)
786 for shiftedWeights, dcrNImage, dcrWeight
in zip(weightsGenerator, dcrNImages, dcrWeights):
787 dcrNImage.array += np.rint(shiftedWeights).astype(dcrNImage.array.dtype)
788 dcrWeight.array += shiftedWeights
790 weightsThreshold = 1.
791 goodPix = dcrWeights[0].array > weightsThreshold
792 for weights
in dcrWeights[1:]:
793 goodPix = (weights.array > weightsThreshold) & goodPix
794 for subfilter
in range(self.config.dcrNumSubfilters):
795 dcrWeights[subfilter].array[goodPix] = 1./dcrWeights[subfilter].array[goodPix]
796 dcrWeights[subfilter].array[~goodPix] = 0.
797 dcrNImages[subfilter].array[~goodPix] = 0
798 return (dcrNImages, dcrWeights)
801 statsCtrl, convergenceMetric,
802 gain, modelWeights, refImage, dcrWeights):
803 """Assemble the DCR coadd for a sub-region.
805 Build a DCR-matched template for each input exposure, then shift the
806 residuals according to the DCR in each subfilter.
807 Stack the shifted residuals and apply them as a correction to the
808 solution from the previous iteration.
809 Restrict the new model solutions from varying by more than a factor of
810 `modelClampFactor` from the last solution, and additionally restrict the
811 individual subfilter models from varying by more than a factor of
812 `frequencyClampFactor` from their average.
813 Finally, mitigate potentially oscillating solutions by averaging the new
814 solution with the solution from the previous iteration, weighted by
815 their convergence metric.
819 dcrModels : `lsst.pipe.tasks.DcrModel`
820 Best fit model of the true sky after correcting chromatic effects.
821 subExposures : `dict` of `lsst.afw.image.ExposureF`
822 The pre-loaded exposures for the current subregion.
823 bbox : `lsst.geom.box.Box2I`
824 Bounding box of the subregion to coadd.
825 dcrBBox : `lsst.geom.box.Box2I`
826 Sub-region of the coadd which includes a buffer to allow for DCR.
827 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
828 `lsst.daf.persistence.ButlerDataRef`
829 The data references to the input warped exposures.
830 statsCtrl : `lsst.afw.math.StatisticsControl`
831 Statistics control object for coadd
832 convergenceMetric : `float`
833 Quality of fit metric for the matched templates of the input images.
834 gain : `float`, optional
835 Relative weight to give the new solution when updating the model.
836 modelWeights : `numpy.ndarray` or `float`
837 A 2D array of weight values that tapers smoothly to zero away from detected sources.
838 Set to a placeholder value of 1.0 if ``self.config.useModelWeights`` is False.
839 refImage : `lsst.afw.image.Image`
840 A reference image used to supply the default pixel values.
841 dcrWeights : `list` of `lsst.afw.image.Image`
842 Per-pixel weights for each subfilter.
843 Equal to 1/(number of unmasked images contributing to each pixel).
845 residualGeneratorList = []
847 for warpExpRef
in warpRefList:
848 visit = warpExpRef.dataId[
"visit"]
849 exposure = subExposures[visit]
850 visitInfo = exposure.getInfo().getVisitInfo()
851 wcs = exposure.getInfo().getWcs()
852 templateImage = dcrModels.buildMatchedTemplate(exposure=exposure,
853 order=self.config.imageInterpOrder,
854 splitSubfilters=self.config.splitSubfilters,
855 splitThreshold=self.config.splitThreshold,
856 amplifyModel=self.config.accelerateModel)
857 residual = exposure.image.array - templateImage.array
859 residual *= exposure.variance.array
863 residualGeneratorList.append(self.dcrResiduals(residual, visitInfo, wcs,
864 dcrModels.effectiveWavelength,
865 dcrModels.bandwidth))
867 dcrSubModelOut = self.newModelFromResidual(dcrModels, residualGeneratorList, dcrBBox, statsCtrl,
869 modelWeights=modelWeights,
871 dcrWeights=dcrWeights)
872 dcrModels.assign(dcrSubModelOut, bbox)
874 def dcrResiduals(self, residual, visitInfo, wcs, effectiveWavelength, bandwidth):
875 """Prepare a residual image for stacking in each subfilter by applying the reverse DCR shifts.
879 residual : `numpy.ndarray`
880 The residual masked image for one exposure,
881 after subtracting the matched template
882 visitInfo : `lsst.afw.image.VisitInfo`
883 Metadata for the exposure.
884 wcs : `lsst.afw.geom.SkyWcs`
885 Coordinate system definition (wcs) for the exposure.
889 residualImage : `numpy.ndarray`
890 The residual image for the next subfilter, shifted for DCR.
894 filteredResidual = ndimage.spline_filter(residual, order=self.config.imageInterpOrder)
897 dcrShift =
calculateDcr(visitInfo, wcs, effectiveWavelength, bandwidth, self.config.dcrNumSubfilters,
898 splitSubfilters=
False)
900 yield applyDcr(filteredResidual, dcr, useInverse=
True, splitSubfilters=
False,
901 doPrefilter=
False, order=self.config.imageInterpOrder)
904 gain, modelWeights, refImage, dcrWeights):
905 """Calculate a new DcrModel from a set of image residuals.
909 dcrModels : `lsst.pipe.tasks.DcrModel`
910 Current model of the true sky after correcting chromatic effects.
911 residualGeneratorList : `generator` of `numpy.ndarray`
912 The residual image for the next subfilter, shifted for DCR.
913 dcrBBox : `lsst.geom.box.Box2I`
914 Sub-region of the coadd which includes a buffer to allow for DCR.
915 statsCtrl : `lsst.afw.math.StatisticsControl`
916 Statistics control object for coadd
918 Relative weight to give the new solution when updating the model.
919 modelWeights : `numpy.ndarray` or `float`
920 A 2D array of weight values that tapers smoothly to zero away from detected sources.
921 Set to a placeholder value of 1.0 if ``self.config.useModelWeights`` is False.
922 refImage : `lsst.afw.image.Image`
923 A reference image used to supply the default pixel values.
924 dcrWeights : `list` of `lsst.afw.image.Image`
925 Per-pixel weights for each subfilter.
926 Equal to 1/(number of unmasked images contributing to each pixel).
930 dcrModel : `lsst.pipe.tasks.DcrModel`
931 New model of the true sky after correcting chromatic effects.
934 for subfilter, model
in enumerate(dcrModels):
935 residualsList = [
next(residualGenerator)
for residualGenerator
in residualGeneratorList]
936 residual = np.sum(residualsList, axis=0)
937 residual *= dcrWeights[subfilter][dcrBBox].array
939 newModel = model[dcrBBox].
clone()
940 newModel.array += residual
942 badPixels = ~np.isfinite(newModel.array)
943 newModel.array[badPixels] = model[dcrBBox].array[badPixels]
944 if self.config.regularizeModelIterations > 0:
945 dcrModels.regularizeModelIter(subfilter, newModel, dcrBBox,
946 self.config.regularizeModelIterations,
947 self.config.regularizationWidth)
948 newModelImages.append(newModel)
949 if self.config.regularizeModelFrequency > 0:
950 dcrModels.regularizeModelFreq(newModelImages, dcrBBox, statsCtrl,
951 self.config.regularizeModelFrequency,
952 self.config.regularizationWidth)
953 dcrModels.conditionDcrModel(newModelImages, dcrBBox, gain=gain)
954 self.applyModelWeights(newModelImages, refImage[dcrBBox], modelWeights)
955 return DcrModel(newModelImages, dcrModels.filter, dcrModels.effectiveWavelength,
956 dcrModels.bandwidth, dcrModels.psf,
957 dcrModels.mask, dcrModels.variance)
960 """Calculate a quality of fit metric for the matched templates.
964 dcrModels : `lsst.pipe.tasks.DcrModel`
965 Best fit model of the true sky after correcting chromatic effects.
966 subExposures : `dict` of `lsst.afw.image.ExposureF`
967 The pre-loaded exposures for the current subregion.
968 bbox : `lsst.geom.box.Box2I`
970 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
971 `lsst.daf.persistence.ButlerDataRef`
972 The data references to the input warped exposures.
973 weightList : `list` of `float`
974 The weight to give each input exposure in the coadd
975 statsCtrl : `lsst.afw.math.StatisticsControl`
976 Statistics control object for coadd
980 convergenceMetric : `float`
981 Quality of fit metric for all input exposures, within the sub-region
983 significanceImage = np.abs(dcrModels.getReferenceImage(bbox))
985 significanceImage += nSigma*dcrModels.calculateNoiseCutoff(dcrModels[1], statsCtrl,
986 bufferSize=self.bufferSize)
987 if np.max(significanceImage) == 0:
988 significanceImage += 1.
992 for warpExpRef, expWeight
in zip(warpRefList, weightList):
993 visit = warpExpRef.dataId[
"visit"]
994 exposure = subExposures[visit][bbox]
995 singleMetric = self.calculateSingleConvergence(dcrModels, exposure, significanceImage, statsCtrl)
996 metric += singleMetric
997 metricList[visit] = singleMetric
999 self.log.
info(
"Individual metrics:\n%s", metricList)
1000 return 1.0
if weight == 0.0
else metric/weight
1003 """Calculate a quality of fit metric for a single matched template.
1007 dcrModels : `lsst.pipe.tasks.DcrModel`
1008 Best fit model of the true sky after correcting chromatic effects.
1009 exposure : `lsst.afw.image.ExposureF`
1010 The input warped exposure to evaluate.
1011 significanceImage : `numpy.ndarray`
1012 Array of weights for each pixel corresponding to its significance
1013 for the convergence calculation.
1014 statsCtrl : `lsst.afw.math.StatisticsControl`
1015 Statistics control object for coadd
1019 convergenceMetric : `float`
1020 Quality of fit metric for one exposure, within the sub-region.
1022 convergeMask = exposure.mask.getPlaneBitMask(self.config.convergenceMaskPlanes)
1023 templateImage = dcrModels.buildMatchedTemplate(exposure=exposure,
1024 order=self.config.imageInterpOrder,
1025 splitSubfilters=self.config.splitSubfilters,
1026 splitThreshold=self.config.splitThreshold,
1027 amplifyModel=self.config.accelerateModel)
1028 diffVals = np.abs(exposure.image.array - templateImage.array)*significanceImage
1029 refVals = np.abs(exposure.image.array + templateImage.array)*significanceImage/2.
1031 finitePixels = np.isfinite(diffVals)
1032 goodMaskPixels = (exposure.mask.array & statsCtrl.getAndMask()) == 0
1033 convergeMaskPixels = exposure.mask.array & convergeMask > 0
1034 usePixels = finitePixels & goodMaskPixels & convergeMaskPixels
1035 if np.sum(usePixels) == 0:
1038 diffUse = diffVals[usePixels]
1039 refUse = refVals[usePixels]
1040 metric = np.sum(diffUse/np.median(diffUse))/np.sum(refUse/np.median(diffUse))
1044 """Add a list of sub-band coadds together.
1048 dcrCoadds : `list` of `lsst.afw.image.ExposureF`
1049 A list of coadd exposures, each exposure containing
1050 the model for one subfilter.
1054 coaddExposure : `lsst.afw.image.ExposureF`
1055 A single coadd exposure that is the sum of the sub-bands.
1057 coaddExposure = dcrCoadds[0].
clone()
1058 for coadd
in dcrCoadds[1:]:
1059 coaddExposure.maskedImage += coadd.maskedImage
1060 return coaddExposure
1062 def fillCoadd(self, dcrModels, skyInfo, warpRefList, weightList, calibration=None, coaddInputs=None,
1063 mask=None, variance=None):
1064 """Create a list of coadd exposures from a list of masked images.
1068 dcrModels : `lsst.pipe.tasks.DcrModel`
1069 Best fit model of the true sky after correcting chromatic effects.
1070 skyInfo : `lsst.pipe.base.Struct`
1071 Patch geometry information, from getSkyInfo
1072 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
1073 `lsst.daf.persistence.ButlerDataRef`
1074 The data references to the input warped exposures.
1075 weightList : `list` of `float`
1076 The weight to give each input exposure in the coadd
1077 calibration : `lsst.afw.Image.PhotoCalib`, optional
1078 Scale factor to set the photometric calibration of an exposure.
1079 coaddInputs : `lsst.afw.Image.CoaddInputs`, optional
1080 A record of the observations that are included in the coadd.
1081 mask : `lsst.afw.image.Mask`, optional
1082 Optional mask to override the values in the final coadd.
1083 variance : `lsst.afw.image.Image`, optional
1084 Optional variance plane to override the values in the final coadd.
1088 dcrCoadds : `list` of `lsst.afw.image.ExposureF`
1089 A list of coadd exposures, each exposure containing
1090 the model for one subfilter.
1093 refModel = dcrModels.getReferenceImage()
1094 for model
in dcrModels:
1095 if self.config.accelerateModel > 1:
1096 model.array = (model.array - refModel)*self.config.accelerateModel + refModel
1097 coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs)
1098 if calibration
is not None:
1099 coaddExposure.setPhotoCalib(calibration)
1100 if coaddInputs
is not None:
1101 coaddExposure.getInfo().setCoaddInputs(coaddInputs)
1103 self.assembleMetadata(coaddExposure, warpRefList, weightList)
1105 coaddExposure.setPsf(dcrModels.psf)
1107 maskedImage = afwImage.MaskedImageF(dcrModels.bbox)
1108 maskedImage.image = model
1109 maskedImage.mask = dcrModels.mask
1110 maskedImage.variance = dcrModels.variance
1111 coaddExposure.setMaskedImage(maskedImage[skyInfo.bbox])
1112 coaddExposure.setPhotoCalib(self.scaleZeroPoint.getPhotoCalib())
1113 if mask
is not None:
1114 coaddExposure.setMask(mask)
1115 if variance
is not None:
1116 coaddExposure.setVariance(variance)
1117 dcrCoadds.append(coaddExposure)
1121 """Calculate the gain to use for the current iteration.
1123 After calculating a new DcrModel, each value is averaged with the
1124 value in the corresponding pixel from the previous iteration. This
1125 reduces oscillating solutions that iterative techniques are plagued by,
1126 and speeds convergence. By far the biggest changes to the model
1127 happen in the first couple iterations, so we can also use a more
1128 aggressive gain later when the model is changing slowly.
1132 convergenceList : `list` of `float`
1133 The quality of fit metric from each previous iteration.
1134 gainList : `list` of `float`
1135 The gains used in each previous iteration: appended with the new
1137 Gains are numbers between ``self.config.baseGain`` and 1.
1142 Relative weight to give the new solution when updating the model.
1143 A value of 1.0 gives equal weight to both solutions.
1148 If ``len(convergenceList) != len(gainList)+1``.
1150 nIter = len(convergenceList)
1151 if nIter != len(gainList) + 1:
1152 raise ValueError(
"convergenceList (%d) must be one element longer than gainList (%d)."
1153 % (len(convergenceList), len(gainList)))
1155 if self.config.baseGain
is None:
1158 baseGain = 1./(self.config.dcrNumSubfilters - 1)
1160 baseGain = self.config.baseGain
1162 if self.config.useProgressiveGain
and nIter > 2:
1170 estFinalConv = [((1 + gainList[i])*convergenceList[i + 1] - convergenceList[i])/gainList[i]
1171 for i
in range(nIter - 1)]
1174 estFinalConv = np.array(estFinalConv)
1175 estFinalConv[estFinalConv < 0] = 0
1177 estFinalConv = np.median(estFinalConv[
max(nIter - 5, 0):])
1178 lastGain = gainList[-1]
1179 lastConv = convergenceList[-2]
1180 newConv = convergenceList[-1]
1185 predictedConv = (estFinalConv*lastGain + lastConv)/(1. + lastGain)
1191 delta = (predictedConv - newConv)/((lastConv - estFinalConv)/(1 + lastGain))
1192 newGain = 1 -
abs(delta)
1194 newGain = (newGain + lastGain)/2.
1195 gain =
max(baseGain, newGain)
1198 gainList.append(gain)
1202 """Build an array that smoothly tapers to 0 away from detected sources.
1206 dcrModels : `lsst.pipe.tasks.DcrModel`
1207 Best fit model of the true sky after correcting chromatic effects.
1208 dcrBBox : `lsst.geom.box.Box2I`
1209 Sub-region of the coadd which includes a buffer to allow for DCR.
1213 weights : `numpy.ndarray` or `float`
1214 A 2D array of weight values that tapers smoothly to zero away from detected sources.
1215 Set to a placeholder value of 1.0 if ``self.config.useModelWeights`` is False.
1220 If ``useModelWeights`` is set and ``modelWeightsWidth`` is negative.
1222 if not self.config.useModelWeights:
1224 if self.config.modelWeightsWidth < 0:
1225 raise ValueError(
"modelWeightsWidth must not be negative if useModelWeights is set")
1226 convergeMask = dcrModels.mask.getPlaneBitMask(self.config.convergenceMaskPlanes)
1227 convergeMaskPixels = dcrModels.mask[dcrBBox].array & convergeMask > 0
1228 weights = np.zeros_like(dcrModels[0][dcrBBox].array)
1229 weights[convergeMaskPixels] = 1.
1230 weights = ndimage.filters.gaussian_filter(weights, self.config.modelWeightsWidth)
1231 weights /= np.max(weights)
1235 """Smoothly replace model pixel values with those from a
1236 reference at locations away from detected sources.
1240 modelImages : `list` of `lsst.afw.image.Image`
1241 The new DCR model images from the current iteration.
1242 The values will be modified in place.
1243 refImage : `lsst.afw.image.MaskedImage`
1244 A reference image used to supply the default pixel values.
1245 modelWeights : `numpy.ndarray` or `float`
1246 A 2D array of weight values that tapers smoothly to zero away from detected sources.
1247 Set to a placeholder value of 1.0 if ``self.config.useModelWeights`` is False.
1249 if self.config.useModelWeights:
1250 for model
in modelImages:
1251 model.array *= modelWeights
1252 model.array += refImage.array*(1. - modelWeights)/self.config.dcrNumSubfilters
1255 """Pre-load sub-regions of a list of exposures.
1259 bbox : `lsst.geom.box.Box2I`
1261 statsCtrl : `lsst.afw.math.StatisticsControl`
1262 Statistics control object for coadd
1263 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
1264 `lsst.daf.persistence.ButlerDataRef`
1265 The data references to the input warped exposures.
1266 imageScalerList : `list` of `lsst.pipe.task.ImageScaler`
1267 The image scalars correct for the zero point of the exposures.
1268 spanSetMaskList : `list` of `dict` containing spanSet lists, or None
1269 Each element is dict with keys = mask plane name to add the spans to
1273 subExposures : `dict`
1274 The `dict` keys are the visit IDs,
1275 and the values are `lsst.afw.image.ExposureF`
1276 The pre-loaded exposures for the current subregion.
1277 The variance plane contains weights, and not the variance
1279 tempExpName = self.getTempExpDatasetName(self.warpType)
1280 zipIterables = zip(warpRefList, imageScalerList, spanSetMaskList)
1282 for warpExpRef, imageScaler, altMaskSpans
in zipIterables:
1283 if isinstance(warpExpRef, DeferredDatasetHandle):
1284 exposure = warpExpRef.get(parameters={
'bbox': bbox})
1286 exposure = warpExpRef.get(tempExpName +
"_sub", bbox=bbox)
1287 visit = warpExpRef.dataId[
"visit"]
1288 if altMaskSpans
is not None:
1289 self.applyAltMaskPlanes(exposure.mask, altMaskSpans)
1290 imageScaler.scaleMaskedImage(exposure.maskedImage)
1292 exposure.variance.array[:, :] = 0.
1294 exposure.variance.array[(exposure.mask.array & statsCtrl.getAndMask()) == 0] = 1.
1297 exposure.image.array[(exposure.mask.array & statsCtrl.getAndMask()) > 0] = 0.
1298 subExposures[visit] = exposure
1302 """Compute the PSF of the coadd from the exposures with the best seeing.
1306 templateCoadd : `lsst.afw.image.ExposureF`
1307 The initial coadd exposure before accounting for DCR.
1308 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
1309 `lsst.daf.persistence.ButlerDataRef`
1310 The data references to the input warped exposures.
1314 psf : `lsst.meas.algorithms.CoaddPsf`
1315 The average PSF of the input exposures with the best seeing.
1317 sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
1318 tempExpName = self.getTempExpDatasetName(self.warpType)
1321 ccds = templateCoadd.getInfo().getCoaddInputs().ccds
1322 psfRefSize = templateCoadd.getPsf().computeShape().getDeterminantRadius()*sigma2fwhm
1323 psfSizes = np.zeros(len(ccds))
1324 ccdVisits = np.array(ccds[
"visit"])
1325 for warpExpRef
in warpRefList:
1326 if isinstance(warpExpRef, DeferredDatasetHandle):
1328 psf = warpExpRef.get(component=
"psf")
1331 psf = warpExpRef.get(tempExpName).getPsf()
1332 visit = warpExpRef.dataId[
"visit"]
1333 psfSize = psf.computeShape().getDeterminantRadius()*sigma2fwhm
1334 psfSizes[ccdVisits == visit] = psfSize
1338 sizeThreshold =
min(np.median(psfSizes), psfRefSize)
1339 goodPsfs = psfSizes <= sizeThreshold
1340 psf = measAlg.CoaddPsf(ccds[goodPsfs], templateCoadd.getWcs(),
1341 self.config.coaddPsf.makeControl())
An integer coordinate rectangle.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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
Extent< int, N > ceil(Extent< double, N > const &input) noexcept
Return the component-wise ceil (round towards more positive).
def applyDcr(image, dcr, useInverse=False, splitSubfilters=False, splitThreshold=0., doPrefilter=True, order=3)
def calculateDcr(visitInfo, wcs, effectiveWavelength, bandwidth, dcrNumSubfilters, splitSubfilters=False)
def run(self, coaddExposures, bbox, wcs)
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs)
def makeSkyInfo(skyMap, tractId, patchId)
def loadSubExposures(self, bbox, statsCtrl, warpRefList, imageScalerList, spanSetMaskList)
def fillCoadd(self, dcrModels, skyInfo, warpRefList, weightList, calibration=None, coaddInputs=None, mask=None, variance=None)
def applyModelWeights(self, modelImages, refImage, modelWeights)
def calculateSingleConvergence(self, dcrModels, exposure, significanceImage, statsCtrl)
def stackCoadd(self, dcrCoadds)
def calculateConvergence(self, dcrModels, subExposures, bbox, warpRefList, weightList, statsCtrl)
def dcrAssembleSubregion(self, dcrModels, subExposures, bbox, dcrBBox, warpRefList, statsCtrl, convergenceMetric, gain, modelWeights, refImage, dcrWeights)
def calculateGain(self, convergenceList, gainList)
def calculateModelWeights(self, dcrModels, dcrBBox)
def newModelFromResidual(self, dcrModels, residualGeneratorList, dcrBBox, statsCtrl, gain, modelWeights, refImage, dcrWeights)
def selectCoaddPsf(self, templateCoadd, warpRefList)
def dcrResiduals(self, residual, visitInfo, wcs, effectiveWavelength, bandwidth)
Angle abs(Angle const &a)