Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Allow configuring more statistical options for assembleCoadd #142

Merged
merged 1 commit into from
Aug 9, 2017
Merged

Allow configuring more statistical options for assembleCoadd #142

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
95 changes: 64 additions & 31 deletions python/lsst/pipe/tasks/assembleCoadd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
import lsst.coadd.utils as coaddUtils
import lsst.pipe.base as pipeBase
import lsst.meas.algorithms as measAlg
import lsst.log as log
from .coaddBase import CoaddBaseTask, SelectDataIdContainer
from .interpImage import InterpImageTask
from .matchBackgrounds import MatchBackgroundsTask
Expand All @@ -57,19 +58,24 @@ class AssembleCoaddConfig(CoaddBaseTask.ConfigClass):
length=2,
default=(2000, 2000),
)
statistic = pexConfig.Field(
dtype=str,
doc="Main stacking statistic for aggregating over the epochs.",
default="MEANCLIP",
)
doSigmaClip = pexConfig.Field(
dtype=bool,
doc="Perform sigma clipped outlier rejection? If False then compute a simple mean.",
default=True,
doc="Perform sigma clipped outlier rejection with MEANCLIP statistic? (DEPRECATED)",
default=False,
)
sigmaClip = pexConfig.Field(
dtype=float,
doc="Sigma for outlier rejection; ignored if doSigmaClip false.",
doc="Sigma for outlier rejection; ignored if non-clipping statistic selected.",
default=3.0,
)
clipIter = pexConfig.Field(
dtype=int,
doc="Number of iterations of outlier rejection; ignored if doSigmaClip false.",
doc="Number of iterations of outlier rejection; ignored if non-clipping statistic selected.",
default=2,
)
scaleZeroPoint = pexConfig.ConfigurableField(
Expand Down Expand Up @@ -155,8 +161,19 @@ def validate(self):
if self.makeDirect and self.makePsfMatched:
raise ValueError("Currently, assembleCoadd can only make either Direct or PsfMatched Coadds "
"at a time. Set either makeDirect or makePsfMatched to False")


if self.doSigmaClip and self.statistic != "MEANCLIP":
log.warn('doSigmaClip deprecated. To replicate behavior, setting statistic to "MEANCLIP"')
self.statistic = "MEANCLIP"
if self.doInterp and self.statistic not in ['MEAN', 'MEDIAN', 'MEANCLIP', 'VARIANCE', 'VARIANCECLIP']:
raise ValueError("Must set doInterp=False for statistic=%s, which does not "
"compute and set a non-zero coadd variance estimate." % (self.statistic))

unstackableStats = ['NOTHING', 'ERROR', 'ORMASK']
if not hasattr(afwMath.Property, self.statistic) or self.statistic in unstackableStats:
stackableStats = [str(k) for k in afwMath.Property.__members__.keys()
if str(k) not in unstackableStats]
raise ValueError("statistic %s is not allowed. Please choose one of %s."
% (self.statistic, stackableStats))
# \addtogroup LSST_task_documentation
# \{
# \page AssembleCoaddTask
Expand Down Expand Up @@ -345,8 +362,7 @@ def run(self, dataRef, selectDataList=[]):

coaddExp = self.assemble(skyInfo, inputData.tempExpRefList, inputData.imageScalerList,
inputData.weightList,
inputData.backgroundInfoList if self.config.doMatchBackgrounds else None,
doClip=self.config.doSigmaClip)
inputData.backgroundInfoList if self.config.doMatchBackgrounds else None)
if self.config.doMatchBackgrounds:
self.addBackgroundMatchingMetadata(coaddExp, inputData.tempExpRefList,
inputData.backgroundInfoList)
Expand Down Expand Up @@ -546,7 +562,7 @@ def backgroundMatching(self, inputData, refExpDataRef=None, refImageScaler=None)
imageScalerList=newScaleList, backgroundInfoList=newBackgroundStructList)

def assemble(self, skyInfo, tempExpRefList, imageScalerList, weightList, bgInfoList=None,
altMaskList=None, doClip=False, mask=None):
altMaskList=None, mask=None):
"""!
\anchor AssembleCoaddTask.assemble_

Expand All @@ -555,16 +571,15 @@ def assemble(self, skyInfo, tempExpRefList, imageScalerList, weightList, bgInfoL
Assemble the coadd using the provided list of coaddTempExps. Since the full coadd covers a patch (a
large area), the assembly is performed over small areas on the image at a time in order to
conserve memory usage. Iterate over subregions within the outer bbox of the patch using
\ref assembleSubregion to mean-stack the corresponding subregions from the coaddTempExps (with outlier
rejection if config.doSigmaClip=True). Set the edge bits the the coadd mask based on the weight map.
\ref assembleSubregion to stack the corresponding subregions from the coaddTempExps with the
statistic specified. Set the edge bits the coadd mask based on the weight map.

\param[in] skyInfo: Patch geometry information, from getSkyInfo
\param[in] tempExpRefList: List of data references to Warps (previously called CoaddTempExps)
\param[in] imageScalerList: List of image scalers
\param[in] weightList: List of weights
\param[in] bgInfoList: List of background data from background matching, or None
\param[in] altMaskList: List of alternate masks to use rather than those stored with tempExp, or None
\param[in] doClip: Use clipping when codding?
\param[in] mask: Mask to ignore when coadding
\return coadded exposure
"""
Expand All @@ -584,10 +599,7 @@ def assemble(self, skyInfo, tempExpRefList, imageScalerList, weightList, bgInfoL
bit = afwImage.Mask.getMaskPlane(plane)
statsCtrl.setMaskPropagationThreshold(bit, threshold)

if doClip:
statsFlags = afwMath.MEANCLIP
else:
statsFlags = afwMath.MEAN
statsFlags = afwMath.stringToStatisticsProperty(self.config.statistic)

if bgInfoList is None:
bgInfoList = [None]*len(tempExpRefList)
Expand Down Expand Up @@ -665,9 +677,11 @@ def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList

For each coaddTempExp, check for (and swap in) an alternative mask if one is passed. If background
matching is enabled, add the background and background variance from each coaddTempExp. Remove mask
planes listed in config.removeMaskPlanes, Finally, mean-stack
the actual exposures using \ref afwMath.statisticsStack "statisticsStack" with outlier rejection if
config.doSigmaClip=True. Assign the stacked subregion back to the coadd.
planes listed in config.removeMaskPlanes, Finally, stack the actual exposures using
\ref afwMath.statisticsStack "statisticsStack" with the statistic specified
by statsFlags. Typically, the statsFlag will be one of afwMath.MEAN for a mean-stack or
afwMath.MEANCLIP for outlier rejection using an N-sigma clipped mean where N and iterations
are specified by statsCtrl. Assign the stacked subregion back to the coadd.

\param[in] coaddExposure: The target image for the coadd
\param[in] bbox: Sub-region to coadd
Expand All @@ -676,7 +690,7 @@ def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList
\param[in] weightList: List of weights
\param[in] bgInfoList: List of background data from background matching
\param[in] altMaskList: List of alternate masks to use rather than those stored with tempExp, or None
\param[in] statsFlags: Statistic for coadd
\param[in] statsFlags: afwMath.Property object for statistic for coadd
\param[in] statsCtrl: Statistics control object for coadd
"""
self.log.debug("Computing coadd over %s", bbox)
Expand Down Expand Up @@ -966,6 +980,19 @@ def setDefaults(self):
self.clipDetection.thresholdType = "pixel_stdev"
self.sigmaClip = 1.5
self.clipIter = 3
self.statistic = "MEAN"

def validate(self):
if self.doSigmaClip:
log.warn("Additional Sigma-clipping not allowed in Safe-clipped Coadds. "
"Ignoring doSigmaClip.")
self.doSigmaClip = False
if self.statistic != "MEAN":
raise ValueError("Only MEAN statistic allowed for final stacking in SafeClipAssembleCoadd "
"(%s chosen). Please set statistic to MEAN."
% (self.statistic))
AssembleCoaddTask.ConfigClass.validate(self)


## \addtogroup LSST_task_documentation
## \{
Expand Down Expand Up @@ -1139,14 +1166,12 @@ def assemble(self, skyInfo, tempExpRefList, imageScalerList, weightList, bgModel
afwDet.setMaskFromFootprintList(maskClipBig, bigFootprints, maskClipValue)
maskClip |= maskClipBig

# Assemble coadd from base class, but ignoring CLIPPED pixels (doClip is false)
# Assemble coadd from base class, but ignoring CLIPPED pixels
badMaskPlanes = self.config.badMaskPlanes[:]
badMaskPlanes.append("CLIPPED")
badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes)
coaddExp = AssembleCoaddTask.assemble(self, skyInfo, tempExpRefList, imageScalerList, weightList,
bgModelList, result.tempExpClipList,
doClip=False,
mask=badPixelMask)
bgModelList, result.tempExpClipList, mask=badPixelMask)

# Set the coadd CLIPPED mask from the footprints since currently pixels that are masked
# do not get propagated
Expand All @@ -1170,13 +1195,21 @@ def buildDifferenceImage(self, skyInfo, tempExpRefList, imageScalerList, weightL
@param bgModelList: List of background models from background matching
@return Difference image of unclipped and clipped coadd wrapped in an Exposure
"""
# Build the unclipped coadd
coaddMean = AssembleCoaddTask.assemble(self, skyInfo, tempExpRefList, imageScalerList, weightList,
bgModelList, doClip=False)

# Build the clipped coadd
coaddClip = AssembleCoaddTask.assemble(self, skyInfo, tempExpRefList, imageScalerList, weightList,
bgModelList, doClip=True)
# Clone and upcast self.config because current self.config is frozen
config = AssembleCoaddConfig()
# getattr necessary because subtasks do not survive Config.toDict()
configIntersection = {k: getattr(self.config, k)
for k, v in self.config.toDict().items() if (k in config.keys())}
config.update(**configIntersection)

# statistic MEAN copied from self.config.statistic, but for clarity explicitly assign
config.statistic = 'MEAN'
task = AssembleCoaddTask(config=config)
coaddMean = task.assemble(skyInfo, tempExpRefList, imageScalerList, weightList, bgModelList)

config.statistic = 'MEANCLIP'
task = AssembleCoaddTask(config=config)
coaddClip = task.assemble(skyInfo, tempExpRefList, imageScalerList, weightList, bgModelList)

coaddDiff = coaddMean.getMaskedImage().Factory(coaddMean.getMaskedImage())
coaddDiff -= coaddClip.getMaskedImage()
Expand Down