1 from __future__
import division, absolute_import
31 from .coaddBase
import CoaddBaseTask, SelectDataIdContainer
32 from .interpImage
import InterpImageTask
33 from .matchBackgrounds
import MatchBackgroundsTask
34 from .scaleZeroPoint
import ScaleZeroPointTask
35 from .coaddHelpers
import groupPatchExposures, getGroupDataRef
37 __all__ = [
"AssembleCoaddTask"]
40 subregionSize = pexConfig.ListField(
42 doc =
"Width, height of stack subregion size; " \
43 "make small enough that a full stack of images will fit into memory at once.",
45 default = (2000, 2000),
47 doSigmaClip = pexConfig.Field(
49 doc =
"Perform sigma clipped outlier rejection? If False then compute a simple mean.",
52 sigmaClip = pexConfig.Field(
54 doc =
"Sigma for outlier rejection; ignored if doSigmaClip false.",
57 clipIter = pexConfig.Field(
59 doc =
"Number of iterations of outlier rejection; ignored if doSigmaClip false.",
62 scaleZeroPoint = pexConfig.ConfigurableField(
63 target = ScaleZeroPointTask,
64 doc =
"Task to adjust the photometric zero point of the coadd temp exposures",
66 doInterp = pexConfig.Field(
67 doc =
"Interpolate over NaN pixels? Also extrapolate, if necessary, but the results are ugly.",
71 interpImage = pexConfig.ConfigurableField(
72 target = InterpImageTask,
73 doc =
"Task to interpolate (and extrapolate) over NaN pixels",
75 matchBackgrounds = pexConfig.ConfigurableField(
76 target = MatchBackgroundsTask,
77 doc =
"Task to match backgrounds",
79 maxMatchResidualRatio = pexConfig.Field(
80 doc =
"Maximum ratio of the mean squared error of the background matching model to the variance " \
81 "of the difference in backgrounds",
85 maxMatchResidualRMS = pexConfig.Field(
86 doc =
"Maximum RMS of residuals of the background offset fit in matchBackgrounds.",
90 doWrite = pexConfig.Field(
91 doc =
"Persist coadd?",
95 doMatchBackgrounds = pexConfig.Field(
96 doc =
"Match backgrounds of coadd temp exposures before coadding them? " \
97 "If False, the coadd temp expsosures must already have been background subtracted or matched",
101 autoReference = pexConfig.Field(
102 doc =
"Automatically select the coadd temp exposure to use as a reference for background matching? " \
103 "Ignored if doMatchBackgrounds false. " \
104 "If False you must specify the reference temp exposure as the data Id",
108 maskPropagationThresholds = pexConfig.DictField(
111 doc = (
"Threshold (in fractional weight) of rejection at which we propagate a mask plane to "
112 "the coadd; that is, we set the mask bit on the coadd if the fraction the rejected frames "
113 "would have contributed exceeds this value."),
114 default = {
"SAT": 0.1},
118 CoaddBaseTask.ConfigClass.setDefaults(self)
123 """Assemble a coadd from a set of coaddTempExp
125 ConfigClass = AssembleCoaddConfig
126 _DefaultName =
"assembleCoadd"
129 CoaddBaseTask.__init__(self, *args, **kwargs)
130 self.makeSubtask(
"interpImage")
131 self.makeSubtask(
"matchBackgrounds")
132 self.makeSubtask(
"scaleZeroPoint")
135 def run(self, dataRef, selectDataList=[]):
136 """Assemble a coadd from a set of coaddTempExp
138 The coadd is computed as a mean with optional outlier rejection.
140 AssumeCoaddTask performs coaddition of "coadd temporary exposures" ("coaddTempExps"). Each
141 coaddTempExp is the size of a patch and contains data for one run, visit or (for a non-mosaic camera)
142 exposure. The coaddTempExps to coadd are selected from the provided selectDataList based on their
143 overlap with the patch specified by dataRef.
145 By default, coaddTempExps contain backgrounds and hence require config.doMatchBackgrounds=True.
146 Background-subtracted coaddTempExps can be coadded by setting config.doMatchBackgrounds=False.
148 When background matching is enabled, the task may be configured to automatically select a reference
149 exposure (config.autoReference=True). If this is not done, then we require that the input dataRef
150 provides access to a coaddTempExp (dataset type coaddName + 'Coadd_tempExp') which is used as the
153 @param dataRef: Data reference defining the patch for coaddition and the reference coaddTempExp
154 (if config.autoReference=False). Used to access the following data products:
155 - [in] self.config.coaddName + "Coadd_skyMap"
156 - [in] self.config.coaddName + "Coadd_tempExp" (optionally)
157 - [out] self.config.coaddName + "Coadd"
158 @param selectDataList[in]: List of data references to coaddTempExp. Data to be coadded will be
159 selected from this list based on overlap with the patch defined by dataRef.
161 @return: a pipeBase.Struct with fields:
162 - coaddExposure: coadded exposure
164 skyInfo = self.getSkyInfo(dataRef)
165 calExpRefList = self.selectExposures(dataRef, skyInfo, selectDataList=selectDataList)
166 if len(calExpRefList) == 0:
167 self.log.warn(
"No exposures to coadd")
169 self.log.info(
"Coadding %d exposures" % len(calExpRefList))
173 self.log.info(
"Found %d %s" % (len(inputData.tempExpRefList), self.getTempExpDatasetName()))
174 if len(inputData.tempExpRefList) == 0:
175 self.log.warn(
"No coadd temporary exposures found")
177 if self.config.doMatchBackgrounds:
180 if len(inputData.tempExpRefList) == 0:
181 self.log.warn(
"No valid background models")
184 coaddExp = self.
assemble(skyInfo, inputData.tempExpRefList, inputData.imageScalerList,
185 inputData.weightList,
186 inputData.backgroundInfoList
if self.config.doMatchBackgrounds
else None)
187 if self.config.doMatchBackgrounds:
189 inputData.backgroundInfoList)
191 if self.config.doInterp:
192 self.interpImage.interpolateOnePlane(maskedImage=coaddExp.getMaskedImage(), planeName=
"NO_DATA")
194 varArray = coaddExp.getMaskedImage().getVariance().getArray()
195 varArray[:] = numpy.where(varArray > 0, varArray, numpy.inf)
197 if self.config.doWrite:
198 self.writeCoaddOutput(dataRef, coaddExp)
200 return pipeBase.Struct(coaddExposure=coaddExp)
203 """Generate list of coaddTempExp data references
205 @param patchRef: Data reference for patch
206 @param calExpRefList: List of data references for input calexps
207 @return List of coaddTempExp data references
209 butler = patchRef.getButler()
211 self.getTempExpDatasetName())
212 tempExpRefList = [
getGroupDataRef(butler, self.getTempExpDatasetName(), g, groupData.keys)
for
213 g
in groupData.groups.keys()]
214 return tempExpRefList
217 """Construct an image scaler for the background reference frame
219 If there is no reference frame ('autoReference'), then this is a no-op
222 @param dataRef: Data reference for the background reference frame, or None
223 @return image scaler, or None
225 if self.config.autoReference:
229 dataset = self.getTempExpDatasetName()
230 if not dataRef.datasetExists(dataset):
231 raise RuntimeError(
"Could not find reference exposure %s %s." % (dataset, dataRef.dataId))
233 refExposure = dataRef.get(self.getTempExpDatasetName(), immediate=
True)
234 refImageScaler = self.scaleZeroPoint.computeImageScaler(
235 exposure = refExposure,
238 return refImageScaler
242 """Prepare the input warps for coaddition
244 This involves measuring weights and constructing image scalers
245 for each of the inputs.
247 @param refList: List of data references to tempExp
249 - tempExprefList: List of data references to tempExp
250 - weightList: List of weightings
251 - imageScalerList: List of image scalers
254 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
255 statsCtrl.setNumIter(self.config.clipIter)
256 statsCtrl.setAndMask(self.getBadPixelMask())
257 statsCtrl.setNanSafe(
True)
265 tempExpName = self.getTempExpDatasetName()
266 for tempExpRef
in refList:
267 if not tempExpRef.datasetExists(tempExpName):
268 self.log.warn(
"Could not find %s %s; skipping it" % (tempExpName, tempExpRef.dataId))
271 tempExp = tempExpRef.get(tempExpName, immediate=
True)
272 maskedImage = tempExp.getMaskedImage()
273 imageScaler = self.scaleZeroPoint.computeImageScaler(
275 dataRef = tempExpRef,
278 imageScaler.scaleMaskedImage(maskedImage)
280 self.log.warn(
"Scaling failed for %s (skipping it): %s" % (tempExpRef.dataId, e))
283 afwMath.MEANCLIP, statsCtrl)
284 meanVar, meanVarErr = statObj.getResult(afwMath.MEANCLIP);
285 weight = 1.0 / float(meanVar)
286 if not numpy.isfinite(weight):
287 self.log.warn(
"Non-finite weight for %s: skipping" % (tempExpRef.dataId,))
289 self.log.info(
"Weight of %s %s = %0.3f" % (tempExpName, tempExpRef.dataId, weight))
294 tempExpRefList.append(tempExpRef)
295 weightList.append(weight)
296 imageScalerList.append(imageScaler)
298 return pipeBase.Struct(tempExpRefList=tempExpRefList, weightList=weightList,
299 imageScalerList=imageScalerList)
303 """Perform background matching on the prepared inputs
305 If no reference is provided, the background matcher will select one.
307 This method returns a new inputData Struct that can replace the original.
309 @param inputData: Struct from prepareInputs() with tempExpRefList, weightList, imageScalerList
310 @param refExpDataRef: Data reference for background reference tempExp, or None
311 @param refImageScaler: Image scaler for background reference tempExp, or None
313 - tempExprefList: List of data references to tempExp
314 - weightList: List of weightings
315 - imageScalerList: List of image scalers
316 - backgroundInfoList: result from background matching
319 backgroundInfoList = self.matchBackgrounds.run(
320 expRefList = inputData.tempExpRefList,
321 imageScalerList = inputData.imageScalerList,
322 refExpDataRef = refExpDataRef
if not self.config.autoReference
else None,
323 refImageScaler = refImageScaler,
324 expDatasetType = self.getTempExpDatasetName(),
327 self.log.fatal(
"Cannot match backgrounds: %s" % (e))
328 raise pipeBase.TaskError(
"Background matching failed.")
331 newTempExpRefList = []
332 newBackgroundStructList = []
336 for tempExpRef, bgInfo, scaler, weight
in zip(inputData.tempExpRefList, backgroundInfoList,
337 inputData.imageScalerList, inputData.weightList):
338 if not bgInfo.isReference:
341 if (bgInfo.backgroundModel
is None):
342 self.log.info(
"No background offset model available for %s: skipping"%(
346 varianceRatio = bgInfo.matchedMSE / bgInfo.diffImVar
348 self.log.info(
"MSE/Var ratio not calculable (%s) for %s: skipping" %
349 (e, tempExpRef.dataId,))
351 if not numpy.isfinite(varianceRatio):
352 self.log.info(
"MSE/Var ratio not finite (%.2f / %.2f) for %s: skipping" %
353 (bgInfo.matchedMSE, bgInfo.diffImVar,
356 elif (varianceRatio > self.config.maxMatchResidualRatio):
357 self.log.info(
"Bad fit. MSE/Var ratio %.2f > %.2f for %s: skipping" % (
358 varianceRatio, self.config.maxMatchResidualRatio, tempExpRef.dataId,))
360 elif ( bgInfo.fitRMS > self.config.maxMatchResidualRMS):
361 self.log.info(
"Bad fit. RMS %.2f > %.2f for %s: skipping" % (
362 bgInfo.fitRMS, self.config.maxMatchResidualRMS, tempExpRef.dataId,))
364 newWeightList.append(1 / (1 / weight + bgInfo.fitRMS**2))
365 newTempExpRefList.append(tempExpRef)
366 newBackgroundStructList.append(bgInfo)
367 newScaleList.append(scaler)
369 return pipeBase.Struct(tempExpRefList=newTempExpRefList, weightList=newWeightList,
370 imageScalerList=newScaleList, backgroundInfoList=newBackgroundStructList)
372 def assemble(self, skyInfo, tempExpRefList, imageScalerList, weightList, bgInfoList=None):
373 """Assemble a coadd from input warps
375 The assembly is performed over small areas on the image at a time, to
376 conserve memory usage.
378 @param skyInfo: Patch geometry information, from getSkyInfo
379 @param tempExpRefList: List of data references to tempExp
380 @param imageScalerList: List of image scalers
381 @param weightList: List of weights
382 @param bgInfoList: List of background data from background matching
383 @return coadded exposure
385 tempExpName = self.getTempExpDatasetName()
386 self.log.info(
"Assembling %s %s" % (len(tempExpRefList), tempExpName))
389 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
390 statsCtrl.setNumIter(self.config.clipIter)
391 statsCtrl.setAndMask(self.getBadPixelMask())
392 statsCtrl.setNanSafe(
True)
393 statsCtrl.setWeighted(
True)
394 statsCtrl.setCalcErrorFromInputVariance(
True)
395 for plane, threshold
in self.config.maskPropagationThresholds.items():
396 bit = afwImage.MaskU.getMaskPlane(plane)
397 statsCtrl.setMaskPropagationThreshold(bit, threshold)
399 if self.config.doSigmaClip:
400 statsFlags = afwMath.MEANCLIP
402 statsFlags = afwMath.MEAN
404 if bgInfoList
is None:
405 bgInfoList = [
None]*len(tempExpRefList)
407 coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs)
408 coaddExposure.setCalib(self.scaleZeroPoint.getCalib())
409 coaddExposure.getInfo().setCoaddInputs(self.inputRecorder.makeCoaddInputs())
411 coaddMaskedImage = coaddExposure.getMaskedImage()
412 subregionSizeArr = self.config.subregionSize
414 for subBBox
in _subBBoxIter(skyInfo.bbox, subregionSize):
417 weightList, bgInfoList, statsFlags, statsCtrl)
419 self.log.fatal(
"Cannot compute coadd %s: %s" % (subBBox, e,))
426 """Set the metadata for the coadd
428 This basic implementation simply sets the filter from the
431 @param coaddExposure: The target image for the coadd
432 @param tempExpRefList: List of data references to tempExp
433 @param weightList: List of weights
435 tempExpName = self.getTempExpDatasetName()
440 coaddInputs = coaddExposure.getInfo().getCoaddInputs()
441 for tempExpRef, weight
in zip(tempExpRefList, weightList):
442 tempExp = tempExpRef.get(tempExpName +
"_sub", bbox=bbox, imageOrigin=
"LOCAL", immediate=
True)
444 coaddExposure.setFilter(tempExp.getFilter())
446 self.inputRecorder.addVisitToCoadd(coaddInputs, tempExp, weight)
447 coaddInputs.visits.sort()
448 if self.config.doPsfMatch:
449 psf = self.config.modelPsf.apply()
452 coaddExposure.setPsf(psf)
453 apCorrMap = measAlg.makeCoaddApCorrMap(coaddInputs.ccds, coaddExposure.getBBox(afwImage.PARENT),
454 coaddExposure.getWcs())
455 coaddExposure.getInfo().setApCorrMap(apCorrMap)
457 def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList, weightList,
458 bgInfoList, statsFlags, statsCtrl):
459 """Assemble the coadd for a sub-region
461 @param coaddExposure: The target image for the coadd
462 @param bbox: Sub-region to coadd
463 @param tempExpRefList: List of data reference to tempExp
464 @param imageScalerList: List of image scalers
465 @param weightList: List of weights
466 @param bgInfoList: List of background data from background matching
467 @param statsFlags: Statistic for coadd
468 @param statsCtrl: Statistics control object for coadd
470 self.log.logdebug(
"Computing coadd over %s" % bbox)
471 tempExpName = self.getTempExpDatasetName()
472 coaddMaskedImage = coaddExposure.getMaskedImage()
473 coaddView = afwImage.MaskedImageF(coaddMaskedImage, bbox, afwImage.PARENT,
False)
474 maskedImageList = afwImage.vectorMaskedImageF()
475 for tempExpRef, imageScaler, bgInfo
in zip(tempExpRefList, imageScalerList, bgInfoList):
476 exposure = tempExpRef.get(tempExpName +
"_sub", bbox=bbox)
477 maskedImage = exposure.getMaskedImage()
478 imageScaler.scaleMaskedImage(maskedImage)
480 if self.config.doMatchBackgrounds
and not bgInfo.isReference:
481 backgroundModel = bgInfo.backgroundModel
482 backgroundImage = backgroundModel.getImage()
if \
483 self.matchBackgrounds.config.usePolynomial
else \
484 backgroundModel.getImageF()
485 backgroundImage.setXY0(coaddMaskedImage.getXY0())
486 maskedImage += backgroundImage.Factory(backgroundImage, bbox, afwImage.PARENT,
False)
487 var = maskedImage.getVariance()
488 var += (bgInfo.fitRMS)**2
490 maskedImageList.append(maskedImage)
492 with self.timer(
"stack"):
494 maskedImageList, statsFlags, statsCtrl, weightList)
496 coaddView <<= coaddSubregion
500 """Add metadata from the background matching to the coadd
502 @param coaddExposure: Coadd
503 @param backgroundInfoList: List of background info, results from background matching
505 self.log.info(
"Adding exposure information to metadata")
506 metadata = coaddExposure.getMetadata()
507 metadata.addString(
"CTExp_SDQA1_DESCRIPTION",
508 "Background matching: Ratio of matchedMSE / diffImVar")
509 for ind, (tempExpRef, backgroundInfo)
in enumerate(zip(tempExpRefList, backgroundInfoList)):
510 tempExpStr =
'&'.join(
'%s=%s' % (k,v)
for k,v
in tempExpRef.dataId.items())
511 if backgroundInfo.isReference:
512 metadata.addString(
"ReferenceExp_ID", tempExpStr)
514 metadata.addString(
"CTExp_ID_%d" % (ind), tempExpStr)
515 metadata.addDouble(
"CTExp_SDQA1_%d" % (ind),
516 backgroundInfo.matchedMSE/backgroundInfo.diffImVar)
517 metadata.addDouble(
"CTExp_SDQA2_%d" % (ind),
518 backgroundInfo.fitRMS)
521 """Create an argument parser
523 parser = pipeBase.ArgumentParser(name=cls._DefaultName)
524 parser.add_id_argument(
"--id", cls.ConfigClass().coaddName +
"Coadd_tempExp",
525 help=
"data ID, e.g. --id tract=12345 patch=1,2",
526 ContainerClass=AssembleCoaddDataIdContainer)
527 parser.add_id_argument(
"--selectId",
"calexp", help=
"data ID, e.g. --selectId visit=6789 ccd=0..9",
528 ContainerClass=SelectDataIdContainer)
534 """Iterate over subregions of a bbox
536 @param[in] bbox: bounding box over which to iterate: afwGeom.Box2I
537 @param[in] subregionSize: size of sub-bboxes
539 @return subBBox: next sub-bounding box of size subregionSize or smaller;
540 each subBBox is contained within bbox, so it may be smaller than subregionSize at the edges of bbox,
541 but it will never be empty
544 raise RuntimeError(
"bbox %s is empty" % (bbox,))
545 if subregionSize[0] < 1
or subregionSize[1] < 1:
546 raise RuntimeError(
"subregionSize %s must be nonzero" % (subregionSize,))
548 for rowShift
in range(0, bbox.getHeight(), subregionSize[1]):
549 for colShift
in range(0, bbox.getWidth(), subregionSize[0]):
552 if subBBox.isEmpty():
553 raise RuntimeError(
"Bug: empty bbox! bbox=%s, subregionSize=%s, colShift=%s, rowShift=%s" % \
554 (bbox, subregionSize, colShift, rowShift))
560 """A version of lsst.pipe.base.DataIdContainer specialized for assembleCoadd.
563 """Make self.refList from self.idList.
565 Interpret the config.doMatchBackgrounds, config.autoReference,
566 and whether a visit/run supplied.
567 If a visit/run is supplied, config.autoReference is automatically set to False.
568 if config.doMatchBackgrounds == false, then a visit/run will be ignored if accidentally supplied.
571 keysCoadd = namespace.butler.getKeys(datasetType=namespace.config.coaddName +
"Coadd",
573 keysCoaddTempExp = namespace.butler.getKeys(datasetType=namespace.config.coaddName +
"Coadd_tempExp",
576 if namespace.config.doMatchBackgrounds:
577 if namespace.config.autoReference:
578 datasetType = namespace.config.coaddName +
"Coadd"
579 validKeys = keysCoadd
581 datasetType = namespace.config.coaddName +
"Coadd_tempExp"
582 validKeys = keysCoaddTempExp
584 datasetType = namespace.config.coaddName +
"Coadd"
585 validKeys = keysCoadd
587 for dataId
in self.idList:
589 for key
in validKeys:
590 if key
not in dataId:
591 raise RuntimeError(
"--id must include " + key)
594 if (key
not in keysCoadd)
and (key
in keysCoaddTempExp):
595 if namespace.config.autoReference:
597 namespace.config.autoReference =
False
598 datasetType = namespace.config.coaddName +
"Coadd_tempExp"
599 print "Switching config.autoReference to False; applies only to background Matching."
602 dataRef = namespace.butler.dataRef(
603 datasetType = datasetType,
606 self.refList.append(dataRef)
def getBackgroundReferenceScaler
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
def addBackgroundMatchingMetadata
An integer coordinate rectangle.
Pass parameters to a Statistics objectA class to pass parameters which control how the stats are calc...
CoaddPsf is the Psf derived to be used for non-PSF-matched Coadd images.
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
lsst::afw::image::MaskedImage< PixelT >::Ptr statisticsStack(lsst::afw::image::MaskedImage< PixelT > const &image, Property flags, char dimension, StatisticsControl const &sctrl=StatisticsControl())
A function to compute statistics on the rows or columns of an image.