LSSTApplications  19.0.0-14-gb0260a2+72efe9b372,20.0.0+7927753e06,20.0.0+8829bf0056,20.0.0+995114c5d2,20.0.0+b6f4b2abd1,20.0.0+bddc4f4cbe,20.0.0-1-g253301a+8829bf0056,20.0.0-1-g2b7511a+0d71a2d77f,20.0.0-1-g5b95a8c+7461dd0434,20.0.0-12-g321c96ea+23efe4bbff,20.0.0-16-gfab17e72e+fdf35455f6,20.0.0-2-g0070d88+ba3ffc8f0b,20.0.0-2-g4dae9ad+ee58a624b3,20.0.0-2-g61b8584+5d3db074ba,20.0.0-2-gb780d76+d529cf1a41,20.0.0-2-ged6426c+226a441f5f,20.0.0-2-gf072044+8829bf0056,20.0.0-2-gf1f7952+ee58a624b3,20.0.0-20-geae50cf+e37fec0aee,20.0.0-25-g3dcad98+544a109665,20.0.0-25-g5eafb0f+ee58a624b3,20.0.0-27-g64178ef+f1f297b00a,20.0.0-3-g4cc78c6+e0676b0dc8,20.0.0-3-g8f21e14+4fd2c12c9a,20.0.0-3-gbd60e8c+187b78b4b8,20.0.0-3-gbecbe05+48431fa087,20.0.0-38-ge4adf513+a12e1f8e37,20.0.0-4-g97dc21a+544a109665,20.0.0-4-gb4befbc+087873070b,20.0.0-4-gf910f65+5d3db074ba,20.0.0-5-gdfe0fee+199202a608,20.0.0-5-gfbfe500+d529cf1a41,20.0.0-6-g64f541c+d529cf1a41,20.0.0-6-g9a5b7a1+a1cd37312e,20.0.0-68-ga3f3dda+5fca18c6a4,20.0.0-9-g4aef684+e18322736b,w.2020.45
LSSTDataManagementBasePackage
assembleCoadd.py
Go to the documentation of this file.
1 # This file is part of pipe_tasks.
2 #
3 # LSST Data Management System
4 # This product includes software developed by the
5 # LSST Project (http://www.lsst.org/).
6 # See COPYRIGHT file at the top of the source tree.
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <https://www.lsstcorp.org/LegalNotices/>.
21 #
22 import os
23 import copy
24 import numpy
25 import warnings
26 import lsst.pex.config as pexConfig
27 import lsst.pex.exceptions as pexExceptions
28 import lsst.geom as geom
29 import lsst.afw.geom as afwGeom
30 import lsst.afw.image as afwImage
31 import lsst.afw.math as afwMath
32 import lsst.afw.table as afwTable
33 import lsst.afw.detection as afwDet
34 import lsst.coadd.utils as coaddUtils
35 import lsst.pipe.base as pipeBase
36 import lsst.meas.algorithms as measAlg
37 import lsst.log as log
38 import lsstDebug
39 import lsst.utils as utils
40 from .coaddBase import CoaddBaseTask, SelectDataIdContainer, makeSkyInfo, makeCoaddSuffix
41 from .interpImage import InterpImageTask
42 from .scaleZeroPoint import ScaleZeroPointTask
43 from .coaddHelpers import groupPatchExposures, getGroupDataRef
44 from .scaleVariance import ScaleVarianceTask
45 from .maskStreaks import MaskStreaksTask
46 from lsst.meas.algorithms import SourceDetectionTask
47 from lsst.daf.butler import DeferredDatasetHandle
48 
49 __all__ = ["AssembleCoaddTask", "AssembleCoaddConnections", "AssembleCoaddConfig",
50  "SafeClipAssembleCoaddTask", "SafeClipAssembleCoaddConfig",
51  "CompareWarpAssembleCoaddTask", "CompareWarpAssembleCoaddConfig"]
52 
53 
54 class AssembleCoaddConnections(pipeBase.PipelineTaskConnections,
55  dimensions=("tract", "patch", "band", "skymap"),
56  defaultTemplates={"inputCoaddName": "deep",
57  "outputCoaddName": "deep",
58  "warpType": "direct",
59  "warpTypeSuffix": "",
60  "fakesType": ""}):
61  inputWarps = pipeBase.connectionTypes.Input(
62  doc=("Input list of warps to be assemebled i.e. stacked."
63  "WarpType (e.g. direct, psfMatched) is controlled by the warpType config parameter"),
64  name="{inputCoaddName}Coadd_{warpType}Warp",
65  storageClass="ExposureF",
66  dimensions=("tract", "patch", "skymap", "visit", "instrument"),
67  deferLoad=True,
68  multiple=True
69  )
70  skyMap = pipeBase.connectionTypes.Input(
71  doc="Input definition of geometry/bbox and projection/wcs for coadded exposures",
72  name="{inputCoaddName}Coadd_skyMap",
73  storageClass="SkyMap",
74  dimensions=("skymap", ),
75  )
76  brightObjectMask = pipeBase.connectionTypes.PrerequisiteInput(
77  doc=("Input Bright Object Mask mask produced with external catalogs to be applied to the mask plane"
78  " BRIGHT_OBJECT."),
79  name="brightObjectMask",
80  storageClass="ObjectMaskCatalog",
81  dimensions=("tract", "patch", "skymap", "band"),
82  )
83  coaddExposure = pipeBase.connectionTypes.Output(
84  doc="Output coadded exposure, produced by stacking input warps",
85  name="{fakesType}{outputCoaddName}Coadd{warpTypeSuffix}",
86  storageClass="ExposureF",
87  dimensions=("tract", "patch", "skymap", "band"),
88  )
89  nImage = pipeBase.connectionTypes.Output(
90  doc="Output image of number of input images per pixel",
91  name="{outputCoaddName}Coadd_nImage",
92  storageClass="ImageU",
93  dimensions=("tract", "patch", "skymap", "band"),
94  )
95 
96  def __init__(self, *, config=None):
97  super().__init__(config=config)
98 
99  # Override the connection's name template with config to replicate Gen2 behavior
100  # This duplicates some of the logic in the base class, due to wanting Gen2 and
101  # Gen3 configs to stay in sync. This should be removed when gen2 is deprecated
102  templateValues = {name: getattr(config.connections, name) for name in self.defaultTemplates}
103  templateValues['warpType'] = config.warpType
104  templateValues['warpTypeSuffix'] = makeCoaddSuffix(config.warpType)
105  if config.hasFakes:
106  templateValues['fakesType'] = "_fakes"
107  self._nameOverrides = {name: getattr(config.connections, name).format(**templateValues)
108  for name in self.allConnections}
109  self._typeNameToVarName = {v: k for k, v in self._nameOverrides.items()}
110  # End code to remove after deprecation
111 
112  if not config.doMaskBrightObjects:
113  self.prerequisiteInputs.remove("brightObjectMask")
114 
115  if not config.doNImage:
116  self.outputs.remove("nImage")
117 
118 
119 class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig,
120  pipelineConnections=AssembleCoaddConnections):
121  """Configuration parameters for the `AssembleCoaddTask`.
122 
123  Notes
124  -----
125  The `doMaskBrightObjects` and `brightObjectMaskName` configuration options
126  only set the bitplane config.brightObjectMaskName. To make this useful you
127  *must* also configure the flags.pixel algorithm, for example by adding
128 
129  .. code-block:: none
130 
131  config.measurement.plugins["base_PixelFlags"].masksFpCenter.append("BRIGHT_OBJECT")
132  config.measurement.plugins["base_PixelFlags"].masksFpAnywhere.append("BRIGHT_OBJECT")
133 
134  to your measureCoaddSources.py and forcedPhotCoadd.py config overrides.
135  """
136  warpType = pexConfig.Field(
137  doc="Warp name: one of 'direct' or 'psfMatched'",
138  dtype=str,
139  default="direct",
140  )
141  subregionSize = pexConfig.ListField(
142  dtype=int,
143  doc="Width, height of stack subregion size; "
144  "make small enough that a full stack of images will fit into memory at once.",
145  length=2,
146  default=(2000, 2000),
147  )
148  statistic = pexConfig.Field(
149  dtype=str,
150  doc="Main stacking statistic for aggregating over the epochs.",
151  default="MEANCLIP",
152  )
153  doSigmaClip = pexConfig.Field(
154  dtype=bool,
155  doc="Perform sigma clipped outlier rejection with MEANCLIP statistic? (DEPRECATED)",
156  default=False,
157  )
158  sigmaClip = pexConfig.Field(
159  dtype=float,
160  doc="Sigma for outlier rejection; ignored if non-clipping statistic selected.",
161  default=3.0,
162  )
163  clipIter = pexConfig.Field(
164  dtype=int,
165  doc="Number of iterations of outlier rejection; ignored if non-clipping statistic selected.",
166  default=2,
167  )
168  calcErrorFromInputVariance = pexConfig.Field(
169  dtype=bool,
170  doc="Calculate coadd variance from input variance by stacking statistic."
171  "Passed to StatisticsControl.setCalcErrorFromInputVariance()",
172  default=True,
173  )
174  scaleZeroPoint = pexConfig.ConfigurableField(
175  target=ScaleZeroPointTask,
176  doc="Task to adjust the photometric zero point of the coadd temp exposures",
177  )
178  doInterp = pexConfig.Field(
179  doc="Interpolate over NaN pixels? Also extrapolate, if necessary, but the results are ugly.",
180  dtype=bool,
181  default=True,
182  )
183  interpImage = pexConfig.ConfigurableField(
184  target=InterpImageTask,
185  doc="Task to interpolate (and extrapolate) over NaN pixels",
186  )
187  doWrite = pexConfig.Field(
188  doc="Persist coadd?",
189  dtype=bool,
190  default=True,
191  )
192  doNImage = pexConfig.Field(
193  doc="Create image of number of contributing exposures for each pixel",
194  dtype=bool,
195  default=False,
196  )
197  doUsePsfMatchedPolygons = pexConfig.Field(
198  doc="Use ValidPolygons from shrunk Psf-Matched Calexps? Should be set to True by CompareWarp only.",
199  dtype=bool,
200  default=False,
201  )
202  maskPropagationThresholds = pexConfig.DictField(
203  keytype=str,
204  itemtype=float,
205  doc=("Threshold (in fractional weight) of rejection at which we propagate a mask plane to "
206  "the coadd; that is, we set the mask bit on the coadd if the fraction the rejected frames "
207  "would have contributed exceeds this value."),
208  default={"SAT": 0.1},
209  )
210  removeMaskPlanes = pexConfig.ListField(dtype=str, default=["NOT_DEBLENDED"],
211  doc="Mask planes to remove before coadding")
212  doMaskBrightObjects = pexConfig.Field(dtype=bool, default=False,
213  doc="Set mask and flag bits for bright objects?")
214  brightObjectMaskName = pexConfig.Field(dtype=str, default="BRIGHT_OBJECT",
215  doc="Name of mask bit used for bright objects")
216  coaddPsf = pexConfig.ConfigField(
217  doc="Configuration for CoaddPsf",
218  dtype=measAlg.CoaddPsfConfig,
219  )
220  doAttachTransmissionCurve = pexConfig.Field(
221  dtype=bool, default=False, optional=False,
222  doc=("Attach a piecewise TransmissionCurve for the coadd? "
223  "(requires all input Exposures to have TransmissionCurves).")
224  )
225  hasFakes = pexConfig.Field(
226  dtype=bool,
227  default=False,
228  doc="Should be set to True if fake sources have been inserted into the input data."
229  )
230 
231  def setDefaults(self):
232  super().setDefaults()
233  self.badMaskPlanes = ["NO_DATA", "BAD", "SAT", "EDGE"]
234 
235  def validate(self):
236  super().validate()
237  if self.doPsfMatch:
238  # Backwards compatibility.
239  # Configs do not have loggers
240  log.warn("Config doPsfMatch deprecated. Setting warpType='psfMatched'")
241  self.warpType = 'psfMatched'
242  if self.doSigmaClip and self.statistic != "MEANCLIP":
243  log.warn('doSigmaClip deprecated. To replicate behavior, setting statistic to "MEANCLIP"')
244  self.statistic = "MEANCLIP"
245  if self.doInterp and self.statistic not in ['MEAN', 'MEDIAN', 'MEANCLIP', 'VARIANCE', 'VARIANCECLIP']:
246  raise ValueError("Must set doInterp=False for statistic=%s, which does not "
247  "compute and set a non-zero coadd variance estimate." % (self.statistic))
248 
249  unstackableStats = ['NOTHING', 'ERROR', 'ORMASK']
250  if not hasattr(afwMath.Property, self.statistic) or self.statistic in unstackableStats:
251  stackableStats = [str(k) for k in afwMath.Property.__members__.keys()
252  if str(k) not in unstackableStats]
253  raise ValueError("statistic %s is not allowed. Please choose one of %s."
254  % (self.statistic, stackableStats))
255 
256 
257 class AssembleCoaddTask(CoaddBaseTask, pipeBase.PipelineTask):
258  """Assemble a coadded image from a set of warps (coadded temporary exposures).
259 
260  We want to assemble a coadded image from a set of Warps (also called
261  coadded temporary exposures or ``coaddTempExps``).
262  Each input Warp covers a patch on the sky and corresponds to a single
263  run/visit/exposure of the covered patch. We provide the task with a list
264  of Warps (``selectDataList``) from which it selects Warps that cover the
265  specified patch (pointed at by ``dataRef``).
266  Each Warp that goes into a coadd will typically have an independent
267  photometric zero-point. Therefore, we must scale each Warp to set it to
268  a common photometric zeropoint. WarpType may be one of 'direct' or
269  'psfMatched', and the boolean configs `config.makeDirect` and
270  `config.makePsfMatched` set which of the warp types will be coadded.
271  The coadd is computed as a mean with optional outlier rejection.
272  Criteria for outlier rejection are set in `AssembleCoaddConfig`.
273  Finally, Warps can have bad 'NaN' pixels which received no input from the
274  source calExps. We interpolate over these bad (NaN) pixels.
275 
276  `AssembleCoaddTask` uses several sub-tasks. These are
277 
278  - `ScaleZeroPointTask`
279  - create and use an ``imageScaler`` object to scale the photometric zeropoint for each Warp
280  - `InterpImageTask`
281  - interpolate across bad pixels (NaN) in the final coadd
282 
283  You can retarget these subtasks if you wish.
284 
285  Notes
286  -----
287  The `lsst.pipe.base.cmdLineTask.CmdLineTask` interface supports a
288  flag ``-d`` to import ``debug.py`` from your ``PYTHONPATH``; see
289  `baseDebug` for more about ``debug.py`` files. `AssembleCoaddTask` has
290  no debug variables of its own. Some of the subtasks may support debug
291  variables. See the documentation for the subtasks for further information.
292 
293  Examples
294  --------
295  `AssembleCoaddTask` assembles a set of warped images into a coadded image.
296  The `AssembleCoaddTask` can be invoked by running ``assembleCoadd.py``
297  with the flag '--legacyCoadd'. Usage of assembleCoadd.py expects two
298  inputs: a data reference to the tract patch and filter to be coadded, and
299  a list of Warps to attempt to coadd. These are specified using ``--id`` and
300  ``--selectId``, respectively:
301 
302  .. code-block:: none
303 
304  --id = [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]
305  --selectId [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]
306 
307  Only the Warps that cover the specified tract and patch will be coadded.
308  A list of the available optional arguments can be obtained by calling
309  ``assembleCoadd.py`` with the ``--help`` command line argument:
310 
311  .. code-block:: none
312 
313  assembleCoadd.py --help
314 
315  To demonstrate usage of the `AssembleCoaddTask` in the larger context of
316  multi-band processing, we will generate the HSC-I & -R band coadds from
317  HSC engineering test data provided in the ``ci_hsc`` package. To begin,
318  assuming that the lsst stack has been already set up, we must set up the
319  obs_subaru and ``ci_hsc`` packages. This defines the environment variable
320  ``$CI_HSC_DIR`` and points at the location of the package. The raw HSC
321  data live in the ``$CI_HSC_DIR/raw directory``. To begin assembling the
322  coadds, we must first
323 
324  - processCcd
325  - process the individual ccds in $CI_HSC_RAW to produce calibrated exposures
326  - makeSkyMap
327  - create a skymap that covers the area of the sky present in the raw exposures
328  - makeCoaddTempExp
329  - warp the individual calibrated exposures to the tangent plane of the coadd
330 
331  We can perform all of these steps by running
332 
333  .. code-block:: none
334 
335  $CI_HSC_DIR scons warp-903986 warp-904014 warp-903990 warp-904010 warp-903988
336 
337  This will produce warped exposures for each visit. To coadd the warped
338  data, we call assembleCoadd.py as follows:
339 
340  .. code-block:: none
341 
342  assembleCoadd.py --legacyCoadd $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-I \
343  --selectId visit=903986 ccd=16 --selectId visit=903986 ccd=22 --selectId visit=903986 ccd=23 \
344  --selectId visit=903986 ccd=100 --selectId visit=904014 ccd=1 --selectId visit=904014 ccd=6 \
345  --selectId visit=904014 ccd=12 --selectId visit=903990 ccd=18 --selectId visit=903990 ccd=25 \
346  --selectId visit=904010 ccd=4 --selectId visit=904010 ccd=10 --selectId visit=904010 ccd=100 \
347  --selectId visit=903988 ccd=16 --selectId visit=903988 ccd=17 --selectId visit=903988 ccd=23 \
348  --selectId visit=903988 ccd=24
349 
350  that will process the HSC-I band data. The results are written in
351  ``$CI_HSC_DIR/DATA/deepCoadd-results/HSC-I``.
352 
353  You may also choose to run:
354 
355  .. code-block:: none
356 
357  scons warp-903334 warp-903336 warp-903338 warp-903342 warp-903344 warp-903346
358  assembleCoadd.py --legacyCoadd $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-R \
359  --selectId visit=903334 ccd=16 --selectId visit=903334 ccd=22 --selectId visit=903334 ccd=23 \
360  --selectId visit=903334 ccd=100 --selectId visit=903336 ccd=17 --selectId visit=903336 ccd=24 \
361  --selectId visit=903338 ccd=18 --selectId visit=903338 ccd=25 --selectId visit=903342 ccd=4 \
362  --selectId visit=903342 ccd=10 --selectId visit=903342 ccd=100 --selectId visit=903344 ccd=0 \
363  --selectId visit=903344 ccd=5 --selectId visit=903344 ccd=11 --selectId visit=903346 ccd=1 \
364  --selectId visit=903346 ccd=6 --selectId visit=903346 ccd=12
365 
366  to generate the coadd for the HSC-R band if you are interested in
367  following multiBand Coadd processing as discussed in `pipeTasks_multiBand`
368  (but note that normally, one would use the `SafeClipAssembleCoaddTask`
369  rather than `AssembleCoaddTask` to make the coadd.
370  """
371  ConfigClass = AssembleCoaddConfig
372  _DefaultName = "assembleCoadd"
373 
374  def __init__(self, *args, **kwargs):
375  # TODO: DM-17415 better way to handle previously allowed passed args e.g.`AssembleCoaddTask(config)`
376  if args:
377  argNames = ["config", "name", "parentTask", "log"]
378  kwargs.update({k: v for k, v in zip(argNames, args)})
379  warnings.warn("AssembleCoadd received positional args, and casting them as kwargs: %s. "
380  "PipelineTask will not take positional args" % argNames, FutureWarning)
381 
382  super().__init__(**kwargs)
383  self.makeSubtask("interpImage")
384  self.makeSubtask("scaleZeroPoint")
385 
386  if self.config.doMaskBrightObjects:
387  mask = afwImage.Mask()
388  try:
389  self.brightObjectBitmask = 1 << mask.addMaskPlane(self.config.brightObjectMaskName)
390  except pexExceptions.LsstCppException:
391  raise RuntimeError("Unable to define mask plane for bright objects; planes used are %s" %
392  mask.getMaskPlaneDict().keys())
393  del mask
394 
395  self.warpType = self.config.warpType
396 
397  @utils.inheritDoc(pipeBase.PipelineTask)
398  def runQuantum(self, butlerQC, inputRefs, outputRefs):
399  # Docstring to be formatted with info from PipelineTask.runQuantum
400  """
401  Notes
402  -----
403  Assemble a coadd from a set of Warps.
404 
405  PipelineTask (Gen3) entry point to Coadd a set of Warps.
406  Analogous to `runDataRef`, it prepares all the data products to be
407  passed to `run`, and processes the results before returning a struct
408  of results to be written out. AssembleCoadd cannot fit all Warps in memory.
409  Therefore, its inputs are accessed subregion by subregion
410  by the Gen3 `DeferredDatasetHandle` that is analagous to the Gen2
411  `lsst.daf.persistence.ButlerDataRef`. Any updates to this method should
412  correspond to an update in `runDataRef` while both entry points
413  are used.
414  """
415  inputData = butlerQC.get(inputRefs)
416 
417  # Construct skyInfo expected by run
418  # Do not remove skyMap from inputData in case makeSupplementaryDataGen3 needs it
419  skyMap = inputData["skyMap"]
420  outputDataId = butlerQC.quantum.dataId
421 
422  inputData['skyInfo'] = makeSkyInfo(skyMap,
423  tractId=outputDataId['tract'],
424  patchId=outputDataId['patch'])
425 
426  # Construct list of input Deferred Datasets
427  # These quack a bit like like Gen2 DataRefs
428  warpRefList = inputData['inputWarps']
429  # Perform same middle steps as `runDataRef` does
430  inputs = self.prepareInputs(warpRefList)
431  self.log.info("Found %d %s", len(inputs.tempExpRefList),
432  self.getTempExpDatasetName(self.warpType))
433  if len(inputs.tempExpRefList) == 0:
434  self.log.warn("No coadd temporary exposures found")
435  return
436 
437  supplementaryData = self.makeSupplementaryDataGen3(butlerQC, inputRefs, outputRefs)
438  retStruct = self.run(inputData['skyInfo'], inputs.tempExpRefList, inputs.imageScalerList,
439  inputs.weightList, supplementaryData=supplementaryData)
440 
441  inputData.setdefault('brightObjectMask', None)
442  self.processResults(retStruct.coaddExposure, inputData['brightObjectMask'], outputDataId)
443 
444  if self.config.doWrite:
445  butlerQC.put(retStruct, outputRefs)
446  return retStruct
447 
448  @pipeBase.timeMethod
449  def runDataRef(self, dataRef, selectDataList=None, warpRefList=None):
450  """Assemble a coadd from a set of Warps.
451 
452  Pipebase.CmdlineTask entry point to Coadd a set of Warps.
453  Compute weights to be applied to each Warp and
454  find scalings to match the photometric zeropoint to a reference Warp.
455  Assemble the Warps using `run`. Interpolate over NaNs and
456  optionally write the coadd to disk. Return the coadded exposure.
457 
458  Parameters
459  ----------
460  dataRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
461  Data reference defining the patch for coaddition and the
462  reference Warp (if ``config.autoReference=False``).
463  Used to access the following data products:
464  - ``self.config.coaddName + "Coadd_skyMap"``
465  - ``self.config.coaddName + "Coadd_ + <warpType> + "Warp"`` (optionally)
466  - ``self.config.coaddName + "Coadd"``
467  selectDataList : `list`
468  List of data references to Calexps. Data to be coadded will be
469  selected from this list based on overlap with the patch defined
470  by dataRef, grouped by visit, and converted to a list of data
471  references to warps.
472  warpRefList : `list`
473  List of data references to Warps to be coadded.
474  Note: `warpRefList` is just the new name for `tempExpRefList`.
475 
476  Returns
477  -------
478  retStruct : `lsst.pipe.base.Struct`
479  Result struct with components:
480 
481  - ``coaddExposure``: coadded exposure (``Exposure``).
482  - ``nImage``: exposure count image (``Image``).
483  """
484  if selectDataList and warpRefList:
485  raise RuntimeError("runDataRef received both a selectDataList and warpRefList, "
486  "and which to use is ambiguous. Please pass only one.")
487 
488  skyInfo = self.getSkyInfo(dataRef)
489  if warpRefList is None:
490  calExpRefList = self.selectExposures(dataRef, skyInfo, selectDataList=selectDataList)
491  if len(calExpRefList) == 0:
492  self.log.warn("No exposures to coadd")
493  return
494  self.log.info("Coadding %d exposures", len(calExpRefList))
495 
496  warpRefList = self.getTempExpRefList(dataRef, calExpRefList)
497 
498  inputData = self.prepareInputs(warpRefList)
499  self.log.info("Found %d %s", len(inputData.tempExpRefList),
500  self.getTempExpDatasetName(self.warpType))
501  if len(inputData.tempExpRefList) == 0:
502  self.log.warn("No coadd temporary exposures found")
503  return
504 
505  supplementaryData = self.makeSupplementaryData(dataRef, warpRefList=inputData.tempExpRefList)
506 
507  retStruct = self.run(skyInfo, inputData.tempExpRefList, inputData.imageScalerList,
508  inputData.weightList, supplementaryData=supplementaryData)
509 
510  brightObjects = self.readBrightObjectMasks(dataRef) if self.config.doMaskBrightObjects else None
511  self.processResults(retStruct.coaddExposure, brightObjectMasks=brightObjects, dataId=dataRef.dataId)
512 
513  if self.config.doWrite:
514  if self.getCoaddDatasetName(self.warpType) == "deepCoadd" and self.config.hasFakes:
515  coaddDatasetName = "fakes_" + self.getCoaddDatasetName(self.warpType)
516  else:
517  coaddDatasetName = self.getCoaddDatasetName(self.warpType)
518  self.log.info("Persisting %s" % coaddDatasetName)
519  dataRef.put(retStruct.coaddExposure, coaddDatasetName)
520  if self.config.doNImage and retStruct.nImage is not None:
521  dataRef.put(retStruct.nImage, self.getCoaddDatasetName(self.warpType) + '_nImage')
522 
523  return retStruct
524 
525  def processResults(self, coaddExposure, brightObjectMasks=None, dataId=None):
526  """Interpolate over missing data and mask bright stars.
527 
528  Parameters
529  ----------
530  coaddExposure : `lsst.afw.image.Exposure`
531  The coadded exposure to process.
532  dataRef : `lsst.daf.persistence.ButlerDataRef`
533  Butler data reference for supplementary data.
534  """
535  if self.config.doInterp:
536  self.interpImage.run(coaddExposure.getMaskedImage(), planeName="NO_DATA")
537  # The variance must be positive; work around for DM-3201.
538  varArray = coaddExposure.variance.array
539  with numpy.errstate(invalid="ignore"):
540  varArray[:] = numpy.where(varArray > 0, varArray, numpy.inf)
541 
542  if self.config.doMaskBrightObjects:
543  self.setBrightObjectMasks(coaddExposure, brightObjectMasks, dataId)
544 
545  def makeSupplementaryData(self, dataRef, selectDataList=None, warpRefList=None):
546  """Make additional inputs to run() specific to subclasses (Gen2)
547 
548  Duplicates interface of `runDataRef` method
549  Available to be implemented by subclasses only if they need the
550  coadd dataRef for performing preliminary processing before
551  assembling the coadd.
552 
553  Parameters
554  ----------
555  dataRef : `lsst.daf.persistence.ButlerDataRef`
556  Butler data reference for supplementary data.
557  selectDataList : `list` (optional)
558  Optional List of data references to Calexps.
559  warpRefList : `list` (optional)
560  Optional List of data references to Warps.
561  """
562  return pipeBase.Struct()
563 
564  def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs):
565  """Make additional inputs to run() specific to subclasses (Gen3)
566 
567  Duplicates interface of `runQuantum` method.
568  Available to be implemented by subclasses only if they need the
569  coadd dataRef for performing preliminary processing before
570  assembling the coadd.
571 
572  Parameters
573  ----------
574  butlerQC : `lsst.pipe.base.ButlerQuantumContext`
575  Gen3 Butler object for fetching additional data products before
576  running the Task specialized for quantum being processed
577  inputRefs : `lsst.pipe.base.InputQuantizedConnection`
578  Attributes are the names of the connections describing input dataset types.
579  Values are DatasetRefs that task consumes for corresponding dataset type.
580  DataIds are guaranteed to match data objects in ``inputData``.
581  outputRefs : `lsst.pipe.base.OutputQuantizedConnection`
582  Attributes are the names of the connections describing output dataset types.
583  Values are DatasetRefs that task is to produce
584  for corresponding dataset type.
585  """
586  return pipeBase.Struct()
587 
588  def getTempExpRefList(self, patchRef, calExpRefList):
589  """Generate list data references corresponding to warped exposures
590  that lie within the patch to be coadded.
591 
592  Parameters
593  ----------
594  patchRef : `dataRef`
595  Data reference for patch.
596  calExpRefList : `list`
597  List of data references for input calexps.
598 
599  Returns
600  -------
601  tempExpRefList : `list`
602  List of Warp/CoaddTempExp data references.
603  """
604  butler = patchRef.getButler()
605  groupData = groupPatchExposures(patchRef, calExpRefList, self.getCoaddDatasetName(self.warpType),
606  self.getTempExpDatasetName(self.warpType))
607  tempExpRefList = [getGroupDataRef(butler, self.getTempExpDatasetName(self.warpType),
608  g, groupData.keys) for
609  g in groupData.groups.keys()]
610  return tempExpRefList
611 
612  def prepareInputs(self, refList):
613  """Prepare the input warps for coaddition by measuring the weight for
614  each warp and the scaling for the photometric zero point.
615 
616  Each Warp has its own photometric zeropoint and background variance.
617  Before coadding these Warps together, compute a scale factor to
618  normalize the photometric zeropoint and compute the weight for each Warp.
619 
620  Parameters
621  ----------
622  refList : `list`
623  List of data references to tempExp
624 
625  Returns
626  -------
627  result : `lsst.pipe.base.Struct`
628  Result struct with components:
629 
630  - ``tempExprefList``: `list` of data references to tempExp.
631  - ``weightList``: `list` of weightings.
632  - ``imageScalerList``: `list` of image scalers.
633  """
634  statsCtrl = afwMath.StatisticsControl()
635  statsCtrl.setNumSigmaClip(self.config.sigmaClip)
636  statsCtrl.setNumIter(self.config.clipIter)
637  statsCtrl.setAndMask(self.getBadPixelMask())
638  statsCtrl.setNanSafe(True)
639  # compute tempExpRefList: a list of tempExpRef that actually exist
640  # and weightList: a list of the weight of the associated coadd tempExp
641  # and imageScalerList: a list of scale factors for the associated coadd tempExp
642  tempExpRefList = []
643  weightList = []
644  imageScalerList = []
645  tempExpName = self.getTempExpDatasetName(self.warpType)
646  for tempExpRef in refList:
647  # Gen3's DeferredDatasetHandles are guaranteed to exist and
648  # therefore have no datasetExists() method
649  if not isinstance(tempExpRef, DeferredDatasetHandle):
650  if not tempExpRef.datasetExists(tempExpName):
651  self.log.warn("Could not find %s %s; skipping it", tempExpName, tempExpRef.dataId)
652  continue
653 
654  tempExp = tempExpRef.get(datasetType=tempExpName, immediate=True)
655  # Ignore any input warp that is empty of data
656  if numpy.isnan(tempExp.image.array).all():
657  continue
658  maskedImage = tempExp.getMaskedImage()
659  imageScaler = self.scaleZeroPoint.computeImageScaler(
660  exposure=tempExp,
661  dataRef=tempExpRef,
662  )
663  try:
664  imageScaler.scaleMaskedImage(maskedImage)
665  except Exception as e:
666  self.log.warn("Scaling failed for %s (skipping it): %s", tempExpRef.dataId, e)
667  continue
668  statObj = afwMath.makeStatistics(maskedImage.getVariance(), maskedImage.getMask(),
669  afwMath.MEANCLIP, statsCtrl)
670  meanVar, meanVarErr = statObj.getResult(afwMath.MEANCLIP)
671  weight = 1.0 / float(meanVar)
672  if not numpy.isfinite(weight):
673  self.log.warn("Non-finite weight for %s: skipping", tempExpRef.dataId)
674  continue
675  self.log.info("Weight of %s %s = %0.3f", tempExpName, tempExpRef.dataId, weight)
676 
677  del maskedImage
678  del tempExp
679 
680  tempExpRefList.append(tempExpRef)
681  weightList.append(weight)
682  imageScalerList.append(imageScaler)
683 
684  return pipeBase.Struct(tempExpRefList=tempExpRefList, weightList=weightList,
685  imageScalerList=imageScalerList)
686 
687  def prepareStats(self, mask=None):
688  """Prepare the statistics for coadding images.
689 
690  Parameters
691  ----------
692  mask : `int`, optional
693  Bit mask value to exclude from coaddition.
694 
695  Returns
696  -------
697  stats : `lsst.pipe.base.Struct`
698  Statistics structure with the following fields:
699 
700  - ``statsCtrl``: Statistics control object for coadd
701  (`lsst.afw.math.StatisticsControl`)
702  - ``statsFlags``: Statistic for coadd (`lsst.afw.math.Property`)
703  """
704  if mask is None:
705  mask = self.getBadPixelMask()
706  statsCtrl = afwMath.StatisticsControl()
707  statsCtrl.setNumSigmaClip(self.config.sigmaClip)
708  statsCtrl.setNumIter(self.config.clipIter)
709  statsCtrl.setAndMask(mask)
710  statsCtrl.setNanSafe(True)
711  statsCtrl.setWeighted(True)
712  statsCtrl.setCalcErrorFromInputVariance(self.config.calcErrorFromInputVariance)
713  for plane, threshold in self.config.maskPropagationThresholds.items():
714  bit = afwImage.Mask.getMaskPlane(plane)
715  statsCtrl.setMaskPropagationThreshold(bit, threshold)
716  statsFlags = afwMath.stringToStatisticsProperty(self.config.statistic)
717  return pipeBase.Struct(ctrl=statsCtrl, flags=statsFlags)
718 
719  @pipeBase.timeMethod
720  def run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
721  altMaskList=None, mask=None, supplementaryData=None):
722  """Assemble a coadd from input warps
723 
724  Assemble the coadd using the provided list of coaddTempExps. Since
725  the full coadd covers a patch (a large area), the assembly is
726  performed over small areas on the image at a time in order to
727  conserve memory usage. Iterate over subregions within the outer
728  bbox of the patch using `assembleSubregion` to stack the corresponding
729  subregions from the coaddTempExps with the statistic specified.
730  Set the edge bits the coadd mask based on the weight map.
731 
732  Parameters
733  ----------
734  skyInfo : `lsst.pipe.base.Struct`
735  Struct with geometric information about the patch.
736  tempExpRefList : `list`
737  List of data references to Warps (previously called CoaddTempExps).
738  imageScalerList : `list`
739  List of image scalers.
740  weightList : `list`
741  List of weights
742  altMaskList : `list`, optional
743  List of alternate masks to use rather than those stored with
744  tempExp.
745  mask : `int`, optional
746  Bit mask value to exclude from coaddition.
747  supplementaryData : lsst.pipe.base.Struct, optional
748  Struct with additional data products needed to assemble coadd.
749  Only used by subclasses that implement `makeSupplementaryData`
750  and override `run`.
751 
752  Returns
753  -------
754  result : `lsst.pipe.base.Struct`
755  Result struct with components:
756 
757  - ``coaddExposure``: coadded exposure (``lsst.afw.image.Exposure``).
758  - ``nImage``: exposure count image (``lsst.afw.image.Image``), if requested.
759  - ``warpRefList``: input list of refs to the warps (
760  ``lsst.daf.butler.DeferredDatasetHandle`` or
761  ``lsst.daf.persistence.ButlerDataRef``)
762  (unmodified)
763  - ``imageScalerList``: input list of image scalers (unmodified)
764  - ``weightList``: input list of weights (unmodified)
765  """
766  tempExpName = self.getTempExpDatasetName(self.warpType)
767  self.log.info("Assembling %s %s", len(tempExpRefList), tempExpName)
768  stats = self.prepareStats(mask=mask)
769 
770  if altMaskList is None:
771  altMaskList = [None]*len(tempExpRefList)
772 
773  coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs)
774  coaddExposure.setPhotoCalib(self.scaleZeroPoint.getPhotoCalib())
775  coaddExposure.getInfo().setCoaddInputs(self.inputRecorder.makeCoaddInputs())
776  self.assembleMetadata(coaddExposure, tempExpRefList, weightList)
777  coaddMaskedImage = coaddExposure.getMaskedImage()
778  subregionSizeArr = self.config.subregionSize
779  subregionSize = geom.Extent2I(subregionSizeArr[0], subregionSizeArr[1])
780  # if nImage is requested, create a zero one which can be passed to assembleSubregion
781  if self.config.doNImage:
782  nImage = afwImage.ImageU(skyInfo.bbox)
783  else:
784  nImage = None
785  for subBBox in self._subBBoxIter(skyInfo.bbox, subregionSize):
786  try:
787  self.assembleSubregion(coaddExposure, subBBox, tempExpRefList, imageScalerList,
788  weightList, altMaskList, stats.flags, stats.ctrl,
789  nImage=nImage)
790  except Exception as e:
791  self.log.fatal("Cannot compute coadd %s: %s", subBBox, e)
792 
793  self.setInexactPsf(coaddMaskedImage.getMask())
794  # Despite the name, the following doesn't really deal with "EDGE" pixels: it identifies
795  # pixels that didn't receive any unmasked inputs (as occurs around the edge of the field).
796  coaddUtils.setCoaddEdgeBits(coaddMaskedImage.getMask(), coaddMaskedImage.getVariance())
797  return pipeBase.Struct(coaddExposure=coaddExposure, nImage=nImage,
798  warpRefList=tempExpRefList, imageScalerList=imageScalerList,
799  weightList=weightList)
800 
801  def assembleMetadata(self, coaddExposure, tempExpRefList, weightList):
802  """Set the metadata for the coadd.
803 
804  This basic implementation sets the filter from the first input.
805 
806  Parameters
807  ----------
808  coaddExposure : `lsst.afw.image.Exposure`
809  The target exposure for the coadd.
810  tempExpRefList : `list`
811  List of data references to tempExp.
812  weightList : `list`
813  List of weights.
814  """
815  assert len(tempExpRefList) == len(weightList), "Length mismatch"
816  tempExpName = self.getTempExpDatasetName(self.warpType)
817  # We load a single pixel of each coaddTempExp, because we just want to get at the metadata
818  # (and we need more than just the PropertySet that contains the header), which is not possible
819  # with the current butler (see #2777).
820  bbox = geom.Box2I(coaddExposure.getBBox().getMin(), geom.Extent2I(1, 1))
821 
822  if isinstance(tempExpRefList[0], DeferredDatasetHandle):
823  # Gen 3 API
824  tempExpList = [tempExpRef.get(parameters={'bbox': bbox}) for tempExpRef in tempExpRefList]
825  else:
826  # Gen 2 API. Delete this when Gen 2 retired
827  tempExpList = [tempExpRef.get(tempExpName + "_sub", bbox=bbox, immediate=True)
828  for tempExpRef in tempExpRefList]
829  numCcds = sum(len(tempExp.getInfo().getCoaddInputs().ccds) for tempExp in tempExpList)
830 
831  coaddExposure.setFilter(tempExpList[0].getFilter())
832  coaddInputs = coaddExposure.getInfo().getCoaddInputs()
833  coaddInputs.ccds.reserve(numCcds)
834  coaddInputs.visits.reserve(len(tempExpList))
835 
836  for tempExp, weight in zip(tempExpList, weightList):
837  self.inputRecorder.addVisitToCoadd(coaddInputs, tempExp, weight)
838 
839  if self.config.doUsePsfMatchedPolygons:
840  self.shrinkValidPolygons(coaddInputs)
841 
842  coaddInputs.visits.sort()
843  if self.warpType == "psfMatched":
844  # The modelPsf BBox for a psfMatchedWarp/coaddTempExp was dynamically defined by
845  # ModelPsfMatchTask as the square box bounding its spatially-variable, pre-matched WarpedPsf.
846  # Likewise, set the PSF of a PSF-Matched Coadd to the modelPsf
847  # having the maximum width (sufficient because square)
848  modelPsfList = [tempExp.getPsf() for tempExp in tempExpList]
849  modelPsfWidthList = [modelPsf.computeBBox().getWidth() for modelPsf in modelPsfList]
850  psf = modelPsfList[modelPsfWidthList.index(max(modelPsfWidthList))]
851  else:
852  psf = measAlg.CoaddPsf(coaddInputs.ccds, coaddExposure.getWcs(),
853  self.config.coaddPsf.makeControl())
854  coaddExposure.setPsf(psf)
855  apCorrMap = measAlg.makeCoaddApCorrMap(coaddInputs.ccds, coaddExposure.getBBox(afwImage.PARENT),
856  coaddExposure.getWcs())
857  coaddExposure.getInfo().setApCorrMap(apCorrMap)
858  if self.config.doAttachTransmissionCurve:
859  transmissionCurve = measAlg.makeCoaddTransmissionCurve(coaddExposure.getWcs(), coaddInputs.ccds)
860  coaddExposure.getInfo().setTransmissionCurve(transmissionCurve)
861 
862  def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList, weightList,
863  altMaskList, statsFlags, statsCtrl, nImage=None):
864  """Assemble the coadd for a sub-region.
865 
866  For each coaddTempExp, check for (and swap in) an alternative mask
867  if one is passed. Remove mask planes listed in
868  `config.removeMaskPlanes`. Finally, stack the actual exposures using
869  `lsst.afw.math.statisticsStack` with the statistic specified by
870  statsFlags. Typically, the statsFlag will be one of lsst.afw.math.MEAN for
871  a mean-stack or `lsst.afw.math.MEANCLIP` for outlier rejection using
872  an N-sigma clipped mean where N and iterations are specified by
873  statsCtrl. Assign the stacked subregion back to the coadd.
874 
875  Parameters
876  ----------
877  coaddExposure : `lsst.afw.image.Exposure`
878  The target exposure for the coadd.
879  bbox : `lsst.geom.Box`
880  Sub-region to coadd.
881  tempExpRefList : `list`
882  List of data reference to tempExp.
883  imageScalerList : `list`
884  List of image scalers.
885  weightList : `list`
886  List of weights.
887  altMaskList : `list`
888  List of alternate masks to use rather than those stored with
889  tempExp, or None. Each element is dict with keys = mask plane
890  name to which to add the spans.
891  statsFlags : `lsst.afw.math.Property`
892  Property object for statistic for coadd.
893  statsCtrl : `lsst.afw.math.StatisticsControl`
894  Statistics control object for coadd.
895  nImage : `lsst.afw.image.ImageU`, optional
896  Keeps track of exposure count for each pixel.
897  """
898  self.log.debug("Computing coadd over %s", bbox)
899  tempExpName = self.getTempExpDatasetName(self.warpType)
900  coaddExposure.mask.addMaskPlane("REJECTED")
901  coaddExposure.mask.addMaskPlane("CLIPPED")
902  coaddExposure.mask.addMaskPlane("SENSOR_EDGE")
903  maskMap = self.setRejectedMaskMapping(statsCtrl)
904  clipped = afwImage.Mask.getPlaneBitMask("CLIPPED")
905  maskedImageList = []
906  if nImage is not None:
907  subNImage = afwImage.ImageU(bbox.getWidth(), bbox.getHeight())
908  for tempExpRef, imageScaler, altMask in zip(tempExpRefList, imageScalerList, altMaskList):
909 
910  if isinstance(tempExpRef, DeferredDatasetHandle):
911  # Gen 3 API
912  exposure = tempExpRef.get(parameters={'bbox': bbox})
913  else:
914  # Gen 2 API. Delete this when Gen 2 retired
915  exposure = tempExpRef.get(tempExpName + "_sub", bbox=bbox)
916 
917  maskedImage = exposure.getMaskedImage()
918  mask = maskedImage.getMask()
919  if altMask is not None:
920  self.applyAltMaskPlanes(mask, altMask)
921  imageScaler.scaleMaskedImage(maskedImage)
922 
923  # Add 1 for each pixel which is not excluded by the exclude mask.
924  # In legacyCoadd, pixels may also be excluded by afwMath.statisticsStack.
925  if nImage is not None:
926  subNImage.getArray()[maskedImage.getMask().getArray() & statsCtrl.getAndMask() == 0] += 1
927  if self.config.removeMaskPlanes:
928  self.removeMaskPlanes(maskedImage)
929  maskedImageList.append(maskedImage)
930 
931  with self.timer("stack"):
932  coaddSubregion = afwMath.statisticsStack(maskedImageList, statsFlags, statsCtrl, weightList,
933  clipped, # also set output to CLIPPED if sigma-clipped
934  maskMap)
935  coaddExposure.maskedImage.assign(coaddSubregion, bbox)
936  if nImage is not None:
937  nImage.assign(subNImage, bbox)
938 
939  def removeMaskPlanes(self, maskedImage):
940  """Unset the mask of an image for mask planes specified in the config.
941 
942  Parameters
943  ----------
944  maskedImage : `lsst.afw.image.MaskedImage`
945  The masked image to be modified.
946  """
947  mask = maskedImage.getMask()
948  for maskPlane in self.config.removeMaskPlanes:
949  try:
950  mask &= ~mask.getPlaneBitMask(maskPlane)
952  self.log.debug("Unable to remove mask plane %s: no mask plane with that name was found.",
953  maskPlane)
954 
955  @staticmethod
956  def setRejectedMaskMapping(statsCtrl):
957  """Map certain mask planes of the warps to new planes for the coadd.
958 
959  If a pixel is rejected due to a mask value other than EDGE, NO_DATA,
960  or CLIPPED, set it to REJECTED on the coadd.
961  If a pixel is rejected due to EDGE, set the coadd pixel to SENSOR_EDGE.
962  If a pixel is rejected due to CLIPPED, set the coadd pixel to CLIPPED.
963 
964  Parameters
965  ----------
966  statsCtrl : `lsst.afw.math.StatisticsControl`
967  Statistics control object for coadd
968 
969  Returns
970  -------
971  maskMap : `list` of `tuple` of `int`
972  A list of mappings of mask planes of the warped exposures to
973  mask planes of the coadd.
974  """
975  edge = afwImage.Mask.getPlaneBitMask("EDGE")
976  noData = afwImage.Mask.getPlaneBitMask("NO_DATA")
977  clipped = afwImage.Mask.getPlaneBitMask("CLIPPED")
978  toReject = statsCtrl.getAndMask() & (~noData) & (~edge) & (~clipped)
979  maskMap = [(toReject, afwImage.Mask.getPlaneBitMask("REJECTED")),
980  (edge, afwImage.Mask.getPlaneBitMask("SENSOR_EDGE")),
981  (clipped, clipped)]
982  return maskMap
983 
984  def applyAltMaskPlanes(self, mask, altMaskSpans):
985  """Apply in place alt mask formatted as SpanSets to a mask.
986 
987  Parameters
988  ----------
989  mask : `lsst.afw.image.Mask`
990  Original mask.
991  altMaskSpans : `dict`
992  SpanSet lists to apply. Each element contains the new mask
993  plane name (e.g. "CLIPPED and/or "NO_DATA") as the key,
994  and list of SpanSets to apply to the mask.
995 
996  Returns
997  -------
998  mask : `lsst.afw.image.Mask`
999  Updated mask.
1000  """
1001  if self.config.doUsePsfMatchedPolygons:
1002  if ("NO_DATA" in altMaskSpans) and ("NO_DATA" in self.config.badMaskPlanes):
1003  # Clear away any other masks outside the validPolygons. These pixels are no longer
1004  # contributing to inexact PSFs, and will still be rejected because of NO_DATA
1005  # self.config.doUsePsfMatchedPolygons should be True only in CompareWarpAssemble
1006  # This mask-clearing step must only occur *before* applying the new masks below
1007  for spanSet in altMaskSpans['NO_DATA']:
1008  spanSet.clippedTo(mask.getBBox()).clearMask(mask, self.getBadPixelMask())
1009 
1010  for plane, spanSetList in altMaskSpans.items():
1011  maskClipValue = mask.addMaskPlane(plane)
1012  for spanSet in spanSetList:
1013  spanSet.clippedTo(mask.getBBox()).setMask(mask, 2**maskClipValue)
1014  return mask
1015 
1016  def shrinkValidPolygons(self, coaddInputs):
1017  """Shrink coaddInputs' ccds' ValidPolygons in place.
1018 
1019  Either modify each ccd's validPolygon in place, or if CoaddInputs
1020  does not have a validPolygon, create one from its bbox.
1021 
1022  Parameters
1023  ----------
1024  coaddInputs : `lsst.afw.image.coaddInputs`
1025  Original mask.
1026 
1027  """
1028  for ccd in coaddInputs.ccds:
1029  polyOrig = ccd.getValidPolygon()
1030  validPolyBBox = polyOrig.getBBox() if polyOrig else ccd.getBBox()
1031  validPolyBBox.grow(-self.config.matchingKernelSize//2)
1032  if polyOrig:
1033  validPolygon = polyOrig.intersectionSingle(validPolyBBox)
1034  else:
1035  validPolygon = afwGeom.polygon.Polygon(geom.Box2D(validPolyBBox))
1036  ccd.setValidPolygon(validPolygon)
1037 
1038  def readBrightObjectMasks(self, dataRef):
1039  """Retrieve the bright object masks.
1040 
1041  Returns None on failure.
1042 
1043  Parameters
1044  ----------
1045  dataRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
1046  A Butler dataRef.
1047 
1048  Returns
1049  -------
1050  result : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
1051  Bright object mask from the Butler object, or None if it cannot
1052  be retrieved.
1053  """
1054  try:
1055  return dataRef.get(datasetType="brightObjectMask", immediate=True)
1056  except Exception as e:
1057  self.log.warn("Unable to read brightObjectMask for %s: %s", dataRef.dataId, e)
1058  return None
1059 
1060  def setBrightObjectMasks(self, exposure, brightObjectMasks, dataId=None):
1061  """Set the bright object masks.
1062 
1063  Parameters
1064  ----------
1065  exposure : `lsst.afw.image.Exposure`
1066  Exposure under consideration.
1067  dataId : `lsst.daf.persistence.dataId`
1068  Data identifier dict for patch.
1069  brightObjectMasks : `lsst.afw.table`
1070  Table of bright objects to mask.
1071  """
1072 
1073  if brightObjectMasks is None:
1074  self.log.warn("Unable to apply bright object mask: none supplied")
1075  return
1076  self.log.info("Applying %d bright object masks to %s", len(brightObjectMasks), dataId)
1077  mask = exposure.getMaskedImage().getMask()
1078  wcs = exposure.getWcs()
1079  plateScale = wcs.getPixelScale().asArcseconds()
1080 
1081  for rec in brightObjectMasks:
1082  center = geom.PointI(wcs.skyToPixel(rec.getCoord()))
1083  if rec["type"] == "box":
1084  assert rec["angle"] == 0.0, ("Angle != 0 for mask object %s" % rec["id"])
1085  width = rec["width"].asArcseconds()/plateScale # convert to pixels
1086  height = rec["height"].asArcseconds()/plateScale # convert to pixels
1087 
1088  halfSize = geom.ExtentI(0.5*width, 0.5*height)
1089  bbox = geom.Box2I(center - halfSize, center + halfSize)
1090 
1091  bbox = geom.BoxI(geom.PointI(int(center[0] - 0.5*width), int(center[1] - 0.5*height)),
1092  geom.PointI(int(center[0] + 0.5*width), int(center[1] + 0.5*height)))
1093  spans = afwGeom.SpanSet(bbox)
1094  elif rec["type"] == "circle":
1095  radius = int(rec["radius"].asArcseconds()/plateScale) # convert to pixels
1096  spans = afwGeom.SpanSet.fromShape(radius, offset=center)
1097  else:
1098  self.log.warn("Unexpected region type %s at %s" % rec["type"], center)
1099  continue
1100  spans.clippedTo(mask.getBBox()).setMask(mask, self.brightObjectBitmask)
1101 
1102  def setInexactPsf(self, mask):
1103  """Set INEXACT_PSF mask plane.
1104 
1105  If any of the input images isn't represented in the coadd (due to
1106  clipped pixels or chip gaps), the `CoaddPsf` will be inexact. Flag
1107  these pixels.
1108 
1109  Parameters
1110  ----------
1111  mask : `lsst.afw.image.Mask`
1112  Coadded exposure's mask, modified in-place.
1113  """
1114  mask.addMaskPlane("INEXACT_PSF")
1115  inexactPsf = mask.getPlaneBitMask("INEXACT_PSF")
1116  sensorEdge = mask.getPlaneBitMask("SENSOR_EDGE") # chip edges (so PSF is discontinuous)
1117  clipped = mask.getPlaneBitMask("CLIPPED") # pixels clipped from coadd
1118  rejected = mask.getPlaneBitMask("REJECTED") # pixels rejected from coadd due to masks
1119  array = mask.getArray()
1120  selected = array & (sensorEdge | clipped | rejected) > 0
1121  array[selected] |= inexactPsf
1122 
1123  @classmethod
1124  def _makeArgumentParser(cls):
1125  """Create an argument parser.
1126  """
1127  parser = pipeBase.ArgumentParser(name=cls._DefaultName)
1128  parser.add_id_argument("--id", cls.ConfigClass().coaddName + "Coadd_"
1129  + cls.ConfigClass().warpType + "Warp",
1130  help="data ID, e.g. --id tract=12345 patch=1,2",
1131  ContainerClass=AssembleCoaddDataIdContainer)
1132  parser.add_id_argument("--selectId", "calexp", help="data ID, e.g. --selectId visit=6789 ccd=0..9",
1133  ContainerClass=SelectDataIdContainer)
1134  return parser
1135 
1136  @staticmethod
1137  def _subBBoxIter(bbox, subregionSize):
1138  """Iterate over subregions of a bbox.
1139 
1140  Parameters
1141  ----------
1142  bbox : `lsst.geom.Box2I`
1143  Bounding box over which to iterate.
1144  subregionSize: `lsst.geom.Extent2I`
1145  Size of sub-bboxes.
1146 
1147  Yields
1148  ------
1149  subBBox : `lsst.geom.Box2I`
1150  Next sub-bounding box of size ``subregionSize`` or smaller; each ``subBBox``
1151  is contained within ``bbox``, so it may be smaller than ``subregionSize`` at
1152  the edges of ``bbox``, but it will never be empty.
1153  """
1154  if bbox.isEmpty():
1155  raise RuntimeError("bbox %s is empty" % (bbox,))
1156  if subregionSize[0] < 1 or subregionSize[1] < 1:
1157  raise RuntimeError("subregionSize %s must be nonzero" % (subregionSize,))
1158 
1159  for rowShift in range(0, bbox.getHeight(), subregionSize[1]):
1160  for colShift in range(0, bbox.getWidth(), subregionSize[0]):
1161  subBBox = geom.Box2I(bbox.getMin() + geom.Extent2I(colShift, rowShift), subregionSize)
1162  subBBox.clip(bbox)
1163  if subBBox.isEmpty():
1164  raise RuntimeError("Bug: empty bbox! bbox=%s, subregionSize=%s, "
1165  "colShift=%s, rowShift=%s" %
1166  (bbox, subregionSize, colShift, rowShift))
1167  yield subBBox
1168 
1169 
1170 class AssembleCoaddDataIdContainer(pipeBase.DataIdContainer):
1171  """A version of `lsst.pipe.base.DataIdContainer` specialized for assembleCoadd.
1172  """
1173 
1174  def makeDataRefList(self, namespace):
1175  """Make self.refList from self.idList.
1176 
1177  Parameters
1178  ----------
1179  namespace
1180  Results of parsing command-line (with ``butler`` and ``log`` elements).
1181  """
1182  datasetType = namespace.config.coaddName + "Coadd"
1183  keysCoadd = namespace.butler.getKeys(datasetType=datasetType, level=self.level)
1184 
1185  for dataId in self.idList:
1186  # tract and patch are required
1187  for key in keysCoadd:
1188  if key not in dataId:
1189  raise RuntimeError("--id must include " + key)
1190 
1191  dataRef = namespace.butler.dataRef(
1192  datasetType=datasetType,
1193  dataId=dataId,
1194  )
1195  self.refList.append(dataRef)
1196 
1197 
1198 def countMaskFromFootprint(mask, footprint, bitmask, ignoreMask):
1199  """Function to count the number of pixels with a specific mask in a
1200  footprint.
1201 
1202  Find the intersection of mask & footprint. Count all pixels in the mask
1203  that are in the intersection that have bitmask set but do not have
1204  ignoreMask set. Return the count.
1205 
1206  Parameters
1207  ----------
1208  mask : `lsst.afw.image.Mask`
1209  Mask to define intersection region by.
1210  footprint : `lsst.afw.detection.Footprint`
1211  Footprint to define the intersection region by.
1212  bitmask
1213  Specific mask that we wish to count the number of occurances of.
1214  ignoreMask
1215  Pixels to not consider.
1216 
1217  Returns
1218  -------
1219  result : `int`
1220  Count of number of pixels in footprint with specified mask.
1221  """
1222  bbox = footprint.getBBox()
1223  bbox.clip(mask.getBBox(afwImage.PARENT))
1224  fp = afwImage.Mask(bbox)
1225  subMask = mask.Factory(mask, bbox, afwImage.PARENT)
1226  footprint.spans.setMask(fp, bitmask)
1227  return numpy.logical_and((subMask.getArray() & fp.getArray()) > 0,
1228  (subMask.getArray() & ignoreMask) == 0).sum()
1229 
1230 
1231 class SafeClipAssembleCoaddConfig(AssembleCoaddConfig, pipelineConnections=AssembleCoaddConnections):
1232  """Configuration parameters for the SafeClipAssembleCoaddTask.
1233  """
1234  assembleMeanCoadd = pexConfig.ConfigurableField(
1235  target=AssembleCoaddTask,
1236  doc="Task to assemble an initial Coadd using the MEAN statistic.",
1237  )
1238  assembleMeanClipCoadd = pexConfig.ConfigurableField(
1239  target=AssembleCoaddTask,
1240  doc="Task to assemble an initial Coadd using the MEANCLIP statistic.",
1241  )
1242  clipDetection = pexConfig.ConfigurableField(
1243  target=SourceDetectionTask,
1244  doc="Detect sources on difference between unclipped and clipped coadd")
1245  minClipFootOverlap = pexConfig.Field(
1246  doc="Minimum fractional overlap of clipped footprint with visit DETECTED to be clipped",
1247  dtype=float,
1248  default=0.6
1249  )
1250  minClipFootOverlapSingle = pexConfig.Field(
1251  doc="Minimum fractional overlap of clipped footprint with visit DETECTED to be "
1252  "clipped when only one visit overlaps",
1253  dtype=float,
1254  default=0.5
1255  )
1256  minClipFootOverlapDouble = pexConfig.Field(
1257  doc="Minimum fractional overlap of clipped footprints with visit DETECTED to be "
1258  "clipped when two visits overlap",
1259  dtype=float,
1260  default=0.45
1261  )
1262  maxClipFootOverlapDouble = pexConfig.Field(
1263  doc="Maximum fractional overlap of clipped footprints with visit DETECTED when "
1264  "considering two visits",
1265  dtype=float,
1266  default=0.15
1267  )
1268  minBigOverlap = pexConfig.Field(
1269  doc="Minimum number of pixels in footprint to use DETECTED mask from the single visits "
1270  "when labeling clipped footprints",
1271  dtype=int,
1272  default=100
1273  )
1274 
1275  def setDefaults(self):
1276  """Set default values for clipDetection.
1277 
1278  Notes
1279  -----
1280  The numeric values for these configuration parameters were
1281  empirically determined, future work may further refine them.
1282  """
1283  AssembleCoaddConfig.setDefaults(self)
1284  self.clipDetection.doTempLocalBackground = False
1285  self.clipDetection.reEstimateBackground = False
1286  self.clipDetection.returnOriginalFootprints = False
1287  self.clipDetection.thresholdPolarity = "both"
1288  self.clipDetection.thresholdValue = 2
1289  self.clipDetection.nSigmaToGrow = 2
1290  self.clipDetection.minPixels = 4
1291  self.clipDetection.isotropicGrow = True
1292  self.clipDetection.thresholdType = "pixel_stdev"
1293  self.sigmaClip = 1.5
1294  self.clipIter = 3
1295  self.statistic = "MEAN"
1296  self.assembleMeanCoadd.statistic = 'MEAN'
1297  self.assembleMeanClipCoadd.statistic = 'MEANCLIP'
1298  self.assembleMeanCoadd.doWrite = False
1299  self.assembleMeanClipCoadd.doWrite = False
1300 
1301  def validate(self):
1302  if self.doSigmaClip:
1303  log.warn("Additional Sigma-clipping not allowed in Safe-clipped Coadds. "
1304  "Ignoring doSigmaClip.")
1305  self.doSigmaClip = False
1306  if self.statistic != "MEAN":
1307  raise ValueError("Only MEAN statistic allowed for final stacking in SafeClipAssembleCoadd "
1308  "(%s chosen). Please set statistic to MEAN."
1309  % (self.statistic))
1310  AssembleCoaddTask.ConfigClass.validate(self)
1311 
1312 
1313 class SafeClipAssembleCoaddTask(AssembleCoaddTask):
1314  """Assemble a coadded image from a set of coadded temporary exposures,
1315  being careful to clip & flag areas with potential artifacts.
1316 
1317  In ``AssembleCoaddTask``, we compute the coadd as an clipped mean (i.e.,
1318  we clip outliers). The problem with doing this is that when computing the
1319  coadd PSF at a given location, individual visit PSFs from visits with
1320  outlier pixels contribute to the coadd PSF and cannot be treated correctly.
1321  In this task, we correct for this behavior by creating a new
1322  ``badMaskPlane`` 'CLIPPED'. We populate this plane on the input
1323  coaddTempExps and the final coadd where
1324 
1325  i. difference imaging suggests that there is an outlier and
1326  ii. this outlier appears on only one or two images.
1327 
1328  Such regions will not contribute to the final coadd. Furthermore, any
1329  routine to determine the coadd PSF can now be cognizant of clipped regions.
1330  Note that the algorithm implemented by this task is preliminary and works
1331  correctly for HSC data. Parameter modifications and or considerable
1332  redesigning of the algorithm is likley required for other surveys.
1333 
1334  ``SafeClipAssembleCoaddTask`` uses a ``SourceDetectionTask``
1335  "clipDetection" subtask and also sub-classes ``AssembleCoaddTask``.
1336  You can retarget the ``SourceDetectionTask`` "clipDetection" subtask
1337  if you wish.
1338 
1339  Notes
1340  -----
1341  The `lsst.pipe.base.cmdLineTask.CmdLineTask` interface supports a
1342  flag ``-d`` to import ``debug.py`` from your ``PYTHONPATH``;
1343  see `baseDebug` for more about ``debug.py`` files.
1344  `SafeClipAssembleCoaddTask` has no debug variables of its own.
1345  The ``SourceDetectionTask`` "clipDetection" subtasks may support debug
1346  variables. See the documetation for `SourceDetectionTask` "clipDetection"
1347  for further information.
1348 
1349  Examples
1350  --------
1351  `SafeClipAssembleCoaddTask` assembles a set of warped ``coaddTempExp``
1352  images into a coadded image. The `SafeClipAssembleCoaddTask` is invoked by
1353  running assembleCoadd.py *without* the flag '--legacyCoadd'.
1354 
1355  Usage of ``assembleCoadd.py`` expects a data reference to the tract patch
1356  and filter to be coadded (specified using
1357  '--id = [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]')
1358  along with a list of coaddTempExps to attempt to coadd (specified using
1359  '--selectId [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]').
1360  Only the coaddTempExps that cover the specified tract and patch will be
1361  coadded. A list of the available optional arguments can be obtained by
1362  calling assembleCoadd.py with the --help command line argument:
1363 
1364  .. code-block:: none
1365 
1366  assembleCoadd.py --help
1367 
1368  To demonstrate usage of the `SafeClipAssembleCoaddTask` in the larger
1369  context of multi-band processing, we will generate the HSC-I & -R band
1370  coadds from HSC engineering test data provided in the ci_hsc package.
1371  To begin, assuming that the lsst stack has been already set up, we must
1372  set up the obs_subaru and ci_hsc packages. This defines the environment
1373  variable $CI_HSC_DIR and points at the location of the package. The raw
1374  HSC data live in the ``$CI_HSC_DIR/raw`` directory. To begin assembling
1375  the coadds, we must first
1376 
1377  - ``processCcd``
1378  process the individual ccds in $CI_HSC_RAW to produce calibrated exposures
1379  - ``makeSkyMap``
1380  create a skymap that covers the area of the sky present in the raw exposures
1381  - ``makeCoaddTempExp``
1382  warp the individual calibrated exposures to the tangent plane of the coadd</DD>
1383 
1384  We can perform all of these steps by running
1385 
1386  .. code-block:: none
1387 
1388  $CI_HSC_DIR scons warp-903986 warp-904014 warp-903990 warp-904010 warp-903988
1389 
1390  This will produce warped coaddTempExps for each visit. To coadd the
1391  warped data, we call ``assembleCoadd.py`` as follows:
1392 
1393  .. code-block:: none
1394 
1395  assembleCoadd.py $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-I \
1396  --selectId visit=903986 ccd=16 --selectId visit=903986 ccd=22 --selectId visit=903986 ccd=23 \
1397  --selectId visit=903986 ccd=100--selectId visit=904014 ccd=1 --selectId visit=904014 ccd=6 \
1398  --selectId visit=904014 ccd=12 --selectId visit=903990 ccd=18 --selectId visit=903990 ccd=25 \
1399  --selectId visit=904010 ccd=4 --selectId visit=904010 ccd=10 --selectId visit=904010 ccd=100 \
1400  --selectId visit=903988 ccd=16 --selectId visit=903988 ccd=17 --selectId visit=903988 ccd=23 \
1401  --selectId visit=903988 ccd=24
1402 
1403  This will process the HSC-I band data. The results are written in
1404  ``$CI_HSC_DIR/DATA/deepCoadd-results/HSC-I``.
1405 
1406  You may also choose to run:
1407 
1408  .. code-block:: none
1409 
1410  scons warp-903334 warp-903336 warp-903338 warp-903342 warp-903344 warp-903346 nnn
1411  assembleCoadd.py $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-R --selectId visit=903334 ccd=16 \
1412  --selectId visit=903334 ccd=22 --selectId visit=903334 ccd=23 --selectId visit=903334 ccd=100 \
1413  --selectId visit=903336 ccd=17 --selectId visit=903336 ccd=24 --selectId visit=903338 ccd=18 \
1414  --selectId visit=903338 ccd=25 --selectId visit=903342 ccd=4 --selectId visit=903342 ccd=10 \
1415  --selectId visit=903342 ccd=100 --selectId visit=903344 ccd=0 --selectId visit=903344 ccd=5 \
1416  --selectId visit=903344 ccd=11 --selectId visit=903346 ccd=1 --selectId visit=903346 ccd=6 \
1417  --selectId visit=903346 ccd=12
1418 
1419  to generate the coadd for the HSC-R band if you are interested in following
1420  multiBand Coadd processing as discussed in ``pipeTasks_multiBand``.
1421  """
1422  ConfigClass = SafeClipAssembleCoaddConfig
1423  _DefaultName = "safeClipAssembleCoadd"
1424 
1425  def __init__(self, *args, **kwargs):
1426  AssembleCoaddTask.__init__(self, *args, **kwargs)
1427  schema = afwTable.SourceTable.makeMinimalSchema()
1428  self.makeSubtask("clipDetection", schema=schema)
1429  self.makeSubtask("assembleMeanClipCoadd")
1430  self.makeSubtask("assembleMeanCoadd")
1431 
1432  @utils.inheritDoc(AssembleCoaddTask)
1433  @pipeBase.timeMethod
1434  def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, *args, **kwargs):
1435  """Assemble the coadd for a region.
1436 
1437  Compute the difference of coadds created with and without outlier
1438  rejection to identify coadd pixels that have outlier values in some
1439  individual visits.
1440  Detect clipped regions on the difference image and mark these regions
1441  on the one or two individual coaddTempExps where they occur if there
1442  is significant overlap between the clipped region and a source. This
1443  leaves us with a set of footprints from the difference image that have
1444  been identified as having occured on just one or two individual visits.
1445  However, these footprints were generated from a difference image. It
1446  is conceivable for a large diffuse source to have become broken up
1447  into multiple footprints acrosss the coadd difference in this process.
1448  Determine the clipped region from all overlapping footprints from the
1449  detected sources in each visit - these are big footprints.
1450  Combine the small and big clipped footprints and mark them on a new
1451  bad mask plane.
1452  Generate the coadd using `AssembleCoaddTask.run` without outlier
1453  removal. Clipped footprints will no longer make it into the coadd
1454  because they are marked in the new bad mask plane.
1455 
1456  Notes
1457  -----
1458  args and kwargs are passed but ignored in order to match the call
1459  signature expected by the parent task.
1460  """
1461  exp = self.buildDifferenceImage(skyInfo, tempExpRefList, imageScalerList, weightList)
1462  mask = exp.getMaskedImage().getMask()
1463  mask.addMaskPlane("CLIPPED")
1464 
1465  result = self.detectClip(exp, tempExpRefList)
1466 
1467  self.log.info('Found %d clipped objects', len(result.clipFootprints))
1468 
1469  maskClipValue = mask.getPlaneBitMask("CLIPPED")
1470  maskDetValue = mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE")
1471  # Append big footprints from individual Warps to result.clipSpans
1472  bigFootprints = self.detectClipBig(result.clipSpans, result.clipFootprints, result.clipIndices,
1473  result.detectionFootprints, maskClipValue, maskDetValue,
1474  exp.getBBox())
1475  # Create mask of the current clipped footprints
1476  maskClip = mask.Factory(mask.getBBox(afwImage.PARENT))
1477  afwDet.setMaskFromFootprintList(maskClip, result.clipFootprints, maskClipValue)
1478 
1479  maskClipBig = maskClip.Factory(mask.getBBox(afwImage.PARENT))
1480  afwDet.setMaskFromFootprintList(maskClipBig, bigFootprints, maskClipValue)
1481  maskClip |= maskClipBig
1482 
1483  # Assemble coadd from base class, but ignoring CLIPPED pixels
1484  badMaskPlanes = self.config.badMaskPlanes[:]
1485  badMaskPlanes.append("CLIPPED")
1486  badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes)
1487  return AssembleCoaddTask.run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
1488  result.clipSpans, mask=badPixelMask)
1489 
1490  def buildDifferenceImage(self, skyInfo, tempExpRefList, imageScalerList, weightList):
1491  """Return an exposure that contains the difference between unclipped
1492  and clipped coadds.
1493 
1494  Generate a difference image between clipped and unclipped coadds.
1495  Compute the difference image by subtracting an outlier-clipped coadd
1496  from an outlier-unclipped coadd. Return the difference image.
1497 
1498  Parameters
1499  ----------
1500  skyInfo : `lsst.pipe.base.Struct`
1501  Patch geometry information, from getSkyInfo
1502  tempExpRefList : `list`
1503  List of data reference to tempExp
1504  imageScalerList : `list`
1505  List of image scalers
1506  weightList : `list`
1507  List of weights
1508 
1509  Returns
1510  -------
1511  exp : `lsst.afw.image.Exposure`
1512  Difference image of unclipped and clipped coadd wrapped in an Exposure
1513  """
1514  coaddMean = self.assembleMeanCoadd.run(skyInfo, tempExpRefList,
1515  imageScalerList, weightList).coaddExposure
1516 
1517  coaddClip = self.assembleMeanClipCoadd.run(skyInfo, tempExpRefList,
1518  imageScalerList, weightList).coaddExposure
1519 
1520  coaddDiff = coaddMean.getMaskedImage().Factory(coaddMean.getMaskedImage())
1521  coaddDiff -= coaddClip.getMaskedImage()
1522  exp = afwImage.ExposureF(coaddDiff)
1523  exp.setPsf(coaddMean.getPsf())
1524  return exp
1525 
1526  def detectClip(self, exp, tempExpRefList):
1527  """Detect clipped regions on an exposure and set the mask on the
1528  individual tempExp masks.
1529 
1530  Detect footprints in the difference image after smoothing the
1531  difference image with a Gaussian kernal. Identify footprints that
1532  overlap with one or two input ``coaddTempExps`` by comparing the
1533  computed overlap fraction to thresholds set in the config. A different
1534  threshold is applied depending on the number of overlapping visits
1535  (restricted to one or two). If the overlap exceeds the thresholds,
1536  the footprint is considered "CLIPPED" and is marked as such on the
1537  coaddTempExp. Return a struct with the clipped footprints, the indices
1538  of the ``coaddTempExps`` that end up overlapping with the clipped
1539  footprints, and a list of new masks for the ``coaddTempExps``.
1540 
1541  Parameters
1542  ----------
1543  exp : `lsst.afw.image.Exposure`
1544  Exposure to run detection on.
1545  tempExpRefList : `list`
1546  List of data reference to tempExp.
1547 
1548  Returns
1549  -------
1550  result : `lsst.pipe.base.Struct`
1551  Result struct with components:
1552 
1553  - ``clipFootprints``: list of clipped footprints.
1554  - ``clipIndices``: indices for each ``clippedFootprint`` in
1555  ``tempExpRefList``.
1556  - ``clipSpans``: List of dictionaries containing spanSet lists
1557  to clip. Each element contains the new maskplane name
1558  ("CLIPPED") as the key and list of ``SpanSets`` as the value.
1559  - ``detectionFootprints``: List of DETECTED/DETECTED_NEGATIVE plane
1560  compressed into footprints.
1561  """
1562  mask = exp.getMaskedImage().getMask()
1563  maskDetValue = mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE")
1564  fpSet = self.clipDetection.detectFootprints(exp, doSmooth=True, clearMask=True)
1565  # Merge positive and negative together footprints together
1566  fpSet.positive.merge(fpSet.negative)
1567  footprints = fpSet.positive
1568  self.log.info('Found %d potential clipped objects', len(footprints.getFootprints()))
1569  ignoreMask = self.getBadPixelMask()
1570 
1571  clipFootprints = []
1572  clipIndices = []
1573  artifactSpanSets = [{'CLIPPED': list()} for _ in tempExpRefList]
1574 
1575  # for use by detectClipBig
1576  visitDetectionFootprints = []
1577 
1578  dims = [len(tempExpRefList), len(footprints.getFootprints())]
1579  overlapDetArr = numpy.zeros(dims, dtype=numpy.uint16)
1580  ignoreArr = numpy.zeros(dims, dtype=numpy.uint16)
1581 
1582  # Loop over masks once and extract/store only relevant overlap metrics and detection footprints
1583  for i, warpRef in enumerate(tempExpRefList):
1584  tmpExpMask = warpRef.get(datasetType=self.getTempExpDatasetName(self.warpType),
1585  immediate=True).getMaskedImage().getMask()
1586  maskVisitDet = tmpExpMask.Factory(tmpExpMask, tmpExpMask.getBBox(afwImage.PARENT),
1587  afwImage.PARENT, True)
1588  maskVisitDet &= maskDetValue
1589  visitFootprints = afwDet.FootprintSet(maskVisitDet, afwDet.Threshold(1))
1590  visitDetectionFootprints.append(visitFootprints)
1591 
1592  for j, footprint in enumerate(footprints.getFootprints()):
1593  ignoreArr[i, j] = countMaskFromFootprint(tmpExpMask, footprint, ignoreMask, 0x0)
1594  overlapDetArr[i, j] = countMaskFromFootprint(tmpExpMask, footprint, maskDetValue, ignoreMask)
1595 
1596  # build a list of clipped spans for each visit
1597  for j, footprint in enumerate(footprints.getFootprints()):
1598  nPixel = footprint.getArea()
1599  overlap = [] # hold the overlap with each visit
1600  indexList = [] # index of visit in global list
1601  for i in range(len(tempExpRefList)):
1602  ignore = ignoreArr[i, j]
1603  overlapDet = overlapDetArr[i, j]
1604  totPixel = nPixel - ignore
1605 
1606  # If we have more bad pixels than detection skip
1607  if ignore > overlapDet or totPixel <= 0.5*nPixel or overlapDet == 0:
1608  continue
1609  overlap.append(overlapDet/float(totPixel))
1610  indexList.append(i)
1611 
1612  overlap = numpy.array(overlap)
1613  if not len(overlap):
1614  continue
1615 
1616  keep = False # Should this footprint be marked as clipped?
1617  keepIndex = [] # Which tempExps does the clipped footprint belong to
1618 
1619  # If footprint only has one overlap use a lower threshold
1620  if len(overlap) == 1:
1621  if overlap[0] > self.config.minClipFootOverlapSingle:
1622  keep = True
1623  keepIndex = [0]
1624  else:
1625  # This is the general case where only visit should be clipped
1626  clipIndex = numpy.where(overlap > self.config.minClipFootOverlap)[0]
1627  if len(clipIndex) == 1:
1628  keep = True
1629  keepIndex = [clipIndex[0]]
1630 
1631  # Test if there are clipped objects that overlap two different visits
1632  clipIndex = numpy.where(overlap > self.config.minClipFootOverlapDouble)[0]
1633  if len(clipIndex) == 2 and len(overlap) > 3:
1634  clipIndexComp = numpy.where(overlap <= self.config.minClipFootOverlapDouble)[0]
1635  if numpy.max(overlap[clipIndexComp]) <= self.config.maxClipFootOverlapDouble:
1636  keep = True
1637  keepIndex = clipIndex
1638 
1639  if not keep:
1640  continue
1641 
1642  for index in keepIndex:
1643  globalIndex = indexList[index]
1644  artifactSpanSets[globalIndex]['CLIPPED'].append(footprint.spans)
1645 
1646  clipIndices.append(numpy.array(indexList)[keepIndex])
1647  clipFootprints.append(footprint)
1648 
1649  return pipeBase.Struct(clipFootprints=clipFootprints, clipIndices=clipIndices,
1650  clipSpans=artifactSpanSets, detectionFootprints=visitDetectionFootprints)
1651 
1652  def detectClipBig(self, clipList, clipFootprints, clipIndices, detectionFootprints,
1653  maskClipValue, maskDetValue, coaddBBox):
1654  """Return individual warp footprints for large artifacts and append
1655  them to ``clipList`` in place.
1656 
1657  Identify big footprints composed of many sources in the coadd
1658  difference that may have originated in a large diffuse source in the
1659  coadd. We do this by indentifying all clipped footprints that overlap
1660  significantly with each source in all the coaddTempExps.
1661 
1662  Parameters
1663  ----------
1664  clipList : `list`
1665  List of alt mask SpanSets with clipping information. Modified.
1666  clipFootprints : `list`
1667  List of clipped footprints.
1668  clipIndices : `list`
1669  List of which entries in tempExpClipList each footprint belongs to.
1670  maskClipValue
1671  Mask value of clipped pixels.
1672  maskDetValue
1673  Mask value of detected pixels.
1674  coaddBBox : `lsst.geom.Box`
1675  BBox of the coadd and warps.
1676 
1677  Returns
1678  -------
1679  bigFootprintsCoadd : `list`
1680  List of big footprints
1681  """
1682  bigFootprintsCoadd = []
1683  ignoreMask = self.getBadPixelMask()
1684  for index, (clippedSpans, visitFootprints) in enumerate(zip(clipList, detectionFootprints)):
1685  maskVisitDet = afwImage.MaskX(coaddBBox, 0x0)
1686  for footprint in visitFootprints.getFootprints():
1687  footprint.spans.setMask(maskVisitDet, maskDetValue)
1688 
1689  # build a mask of clipped footprints that are in this visit
1690  clippedFootprintsVisit = []
1691  for foot, clipIndex in zip(clipFootprints, clipIndices):
1692  if index not in clipIndex:
1693  continue
1694  clippedFootprintsVisit.append(foot)
1695  maskVisitClip = maskVisitDet.Factory(maskVisitDet.getBBox(afwImage.PARENT))
1696  afwDet.setMaskFromFootprintList(maskVisitClip, clippedFootprintsVisit, maskClipValue)
1697 
1698  bigFootprintsVisit = []
1699  for foot in visitFootprints.getFootprints():
1700  if foot.getArea() < self.config.minBigOverlap:
1701  continue
1702  nCount = countMaskFromFootprint(maskVisitClip, foot, maskClipValue, ignoreMask)
1703  if nCount > self.config.minBigOverlap:
1704  bigFootprintsVisit.append(foot)
1705  bigFootprintsCoadd.append(foot)
1706 
1707  for footprint in bigFootprintsVisit:
1708  clippedSpans["CLIPPED"].append(footprint.spans)
1709 
1710  return bigFootprintsCoadd
1711 
1712 
1714  psfMatchedWarps = pipeBase.connectionTypes.Input(
1715  doc=("PSF-Matched Warps are required by CompareWarp regardless of the coadd type requested. "
1716  "Only PSF-Matched Warps make sense for image subtraction. "
1717  "Therefore, they must be an additional declared input."),
1718  name="{inputCoaddName}Coadd_psfMatchedWarp",
1719  storageClass="ExposureF",
1720  dimensions=("tract", "patch", "skymap", "visit"),
1721  deferLoad=True,
1722  multiple=True
1723  )
1724  templateCoadd = pipeBase.connectionTypes.Output(
1725  doc=("Model of the static sky, used to find temporal artifacts. Typically a PSF-Matched, "
1726  "sigma-clipped coadd. Written if and only if assembleStaticSkyModel.doWrite=True"),
1727  name="{fakesType}{outputCoaddName}CoaddPsfMatched",
1728  storageClass="ExposureF",
1729  dimensions=("tract", "patch", "skymap", "band"),
1730  )
1731 
1732  def __init__(self, *, config=None):
1733  super().__init__(config=config)
1734  if not config.assembleStaticSkyModel.doWrite:
1735  self.outputs.remove("templateCoadd")
1736  config.validate()
1737 
1738 
1739 class CompareWarpAssembleCoaddConfig(AssembleCoaddConfig,
1740  pipelineConnections=CompareWarpAssembleCoaddConnections):
1741  assembleStaticSkyModel = pexConfig.ConfigurableField(
1742  target=AssembleCoaddTask,
1743  doc="Task to assemble an artifact-free, PSF-matched Coadd to serve as a"
1744  " naive/first-iteration model of the static sky.",
1745  )
1746  detect = pexConfig.ConfigurableField(
1747  target=SourceDetectionTask,
1748  doc="Detect outlier sources on difference between each psfMatched warp and static sky model"
1749  )
1750  detectTemplate = pexConfig.ConfigurableField(
1751  target=SourceDetectionTask,
1752  doc="Detect sources on static sky model. Only used if doPreserveContainedBySource is True"
1753  )
1754  maskStreaks = pexConfig.ConfigurableField(
1755  target=MaskStreaksTask,
1756  doc="Detect streaks on difference between each psfMatched warp and static sky model. Only used if "
1757  "doFilterMorphological is True. Adds a mask plane to an exposure, with the mask plane name set by"
1758  "streakMaskName"
1759  )
1760  streakMaskName = pexConfig.Field(
1761  dtype=str,
1762  default="STREAK",
1763  doc="Name of mask bit used for streaks"
1764  )
1765  maxNumEpochs = pexConfig.Field(
1766  doc="Charactistic maximum local number of epochs/visits in which an artifact candidate can appear "
1767  "and still be masked. The effective maxNumEpochs is a broken linear function of local "
1768  "number of epochs (N): min(maxFractionEpochsLow*N, maxNumEpochs + maxFractionEpochsHigh*N). "
1769  "For each footprint detected on the image difference between the psfMatched warp and static sky "
1770  "model, if a significant fraction of pixels (defined by spatialThreshold) are residuals in more "
1771  "than the computed effective maxNumEpochs, the artifact candidate is deemed persistant rather "
1772  "than transient and not masked.",
1773  dtype=int,
1774  default=2
1775  )
1776  maxFractionEpochsLow = pexConfig.RangeField(
1777  doc="Fraction of local number of epochs (N) to use as effective maxNumEpochs for low N. "
1778  "Effective maxNumEpochs = "
1779  "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)",
1780  dtype=float,
1781  default=0.4,
1782  min=0., max=1.,
1783  )
1784  maxFractionEpochsHigh = pexConfig.RangeField(
1785  doc="Fraction of local number of epochs (N) to use as effective maxNumEpochs for high N. "
1786  "Effective maxNumEpochs = "
1787  "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)",
1788  dtype=float,
1789  default=0.03,
1790  min=0., max=1.,
1791  )
1792  spatialThreshold = pexConfig.RangeField(
1793  doc="Unitless fraction of pixels defining how much of the outlier region has to meet the "
1794  "temporal criteria. If 0, clip all. If 1, clip none.",
1795  dtype=float,
1796  default=0.5,
1797  min=0., max=1.,
1798  inclusiveMin=True, inclusiveMax=True
1799  )
1800  doScaleWarpVariance = pexConfig.Field(
1801  doc="Rescale Warp variance plane using empirical noise?",
1802  dtype=bool,
1803  default=True,
1804  )
1805  scaleWarpVariance = pexConfig.ConfigurableField(
1806  target=ScaleVarianceTask,
1807  doc="Rescale variance on warps",
1808  )
1809  doPreserveContainedBySource = pexConfig.Field(
1810  doc="Rescue artifacts from clipping that completely lie within a footprint detected"
1811  "on the PsfMatched Template Coadd. Replicates a behavior of SafeClip.",
1812  dtype=bool,
1813  default=True,
1814  )
1815  doPrefilterArtifacts = pexConfig.Field(
1816  doc="Ignore artifact candidates that are mostly covered by the bad pixel mask, "
1817  "because they will be excluded anyway. This prevents them from contributing "
1818  "to the outlier epoch count image and potentially being labeled as persistant."
1819  "'Mostly' is defined by the config 'prefilterArtifactsRatio'.",
1820  dtype=bool,
1821  default=True
1822  )
1823  prefilterArtifactsMaskPlanes = pexConfig.ListField(
1824  doc="Prefilter artifact candidates that are mostly covered by these bad mask planes.",
1825  dtype=str,
1826  default=('NO_DATA', 'BAD', 'SAT', 'SUSPECT'),
1827  )
1828  prefilterArtifactsRatio = pexConfig.Field(
1829  doc="Prefilter artifact candidates with less than this fraction overlapping good pixels",
1830  dtype=float,
1831  default=0.05
1832  )
1833  doFilterMorphological = pexConfig.Field(
1834  doc="Filter artifact candidates based on morphological criteria, i.g. those that appear to "
1835  "be streaks.",
1836  dtype=bool,
1837  default=False
1838  )
1839 
1840  def setDefaults(self):
1841  AssembleCoaddConfig.setDefaults(self)
1842  self.statistic = 'MEAN'
1844 
1845  # Real EDGE removed by psfMatched NO_DATA border half the width of the matching kernel
1846  # CompareWarp applies psfMatched EDGE pixels to directWarps before assembling
1847  if "EDGE" in self.badMaskPlanes:
1848  self.badMaskPlanes.remove('EDGE')
1849  self.removeMaskPlanes.append('EDGE')
1850  self.assembleStaticSkyModel.badMaskPlanes = ["NO_DATA", ]
1851  self.assembleStaticSkyModel.warpType = 'psfMatched'
1852  self.assembleStaticSkyModel.connections.warpType = 'psfMatched'
1853  self.assembleStaticSkyModel.statistic = 'MEANCLIP'
1854  self.assembleStaticSkyModel.sigmaClip = 2.5
1855  self.assembleStaticSkyModel.clipIter = 3
1856  self.assembleStaticSkyModel.calcErrorFromInputVariance = False
1857  self.assembleStaticSkyModel.doWrite = False
1858  self.detect.doTempLocalBackground = False
1859  self.detect.reEstimateBackground = False
1860  self.detect.returnOriginalFootprints = False
1861  self.detect.thresholdPolarity = "both"
1862  self.detect.thresholdValue = 5
1863  self.detect.minPixels = 4
1864  self.detect.isotropicGrow = True
1865  self.detect.thresholdType = "pixel_stdev"
1866  self.detect.nSigmaToGrow = 0.4
1867  # The default nSigmaToGrow for SourceDetectionTask is already 2.4,
1868  # Explicitly restating because ratio with detect.nSigmaToGrow matters
1869  self.detectTemplate.nSigmaToGrow = 2.4
1870  self.detectTemplate.doTempLocalBackground = False
1871  self.detectTemplate.reEstimateBackground = False
1872  self.detectTemplate.returnOriginalFootprints = False
1873 
1874  def validate(self):
1875  super().validate()
1876  if self.assembleStaticSkyModel.doNImage:
1877  raise ValueError("No dataset type exists for a PSF-Matched Template N Image."
1878  "Please set assembleStaticSkyModel.doNImage=False")
1879 
1880  if self.assembleStaticSkyModel.doWrite and (self.warpType == self.assembleStaticSkyModel.warpType):
1881  raise ValueError("warpType (%s) == assembleStaticSkyModel.warpType (%s) and will compete for "
1882  "the same dataset name. Please set assembleStaticSkyModel.doWrite to False "
1883  "or warpType to 'direct'. assembleStaticSkyModel.warpType should ways be "
1884  "'PsfMatched'" % (self.warpType, self.assembleStaticSkyModel.warpType))
1885 
1886 
1887 class CompareWarpAssembleCoaddTask(AssembleCoaddTask):
1888  """Assemble a compareWarp coadded image from a set of warps
1889  by masking artifacts detected by comparing PSF-matched warps.
1890 
1891  In ``AssembleCoaddTask``, we compute the coadd as an clipped mean (i.e.,
1892  we clip outliers). The problem with doing this is that when computing the
1893  coadd PSF at a given location, individual visit PSFs from visits with
1894  outlier pixels contribute to the coadd PSF and cannot be treated correctly.
1895  In this task, we correct for this behavior by creating a new badMaskPlane
1896  'CLIPPED' which marks pixels in the individual warps suspected to contain
1897  an artifact. We populate this plane on the input warps by comparing
1898  PSF-matched warps with a PSF-matched median coadd which serves as a
1899  model of the static sky. Any group of pixels that deviates from the
1900  PSF-matched template coadd by more than config.detect.threshold sigma,
1901  is an artifact candidate. The candidates are then filtered to remove
1902  variable sources and sources that are difficult to subtract such as
1903  bright stars. This filter is configured using the config parameters
1904  ``temporalThreshold`` and ``spatialThreshold``. The temporalThreshold is
1905  the maximum fraction of epochs that the deviation can appear in and still
1906  be considered an artifact. The spatialThreshold is the maximum fraction of
1907  pixels in the footprint of the deviation that appear in other epochs
1908  (where other epochs is defined by the temporalThreshold). If the deviant
1909  region meets this criteria of having a significant percentage of pixels
1910  that deviate in only a few epochs, these pixels have the 'CLIPPED' bit
1911  set in the mask. These regions will not contribute to the final coadd.
1912  Furthermore, any routine to determine the coadd PSF can now be cognizant
1913  of clipped regions. Note that the algorithm implemented by this task is
1914  preliminary and works correctly for HSC data. Parameter modifications and
1915  or considerable redesigning of the algorithm is likley required for other
1916  surveys.
1917 
1918  ``CompareWarpAssembleCoaddTask`` sub-classes
1919  ``AssembleCoaddTask`` and instantiates ``AssembleCoaddTask``
1920  as a subtask to generate the TemplateCoadd (the model of the static sky).
1921 
1922  Notes
1923  -----
1924  The `lsst.pipe.base.cmdLineTask.CmdLineTask` interface supports a
1925  flag ``-d`` to import ``debug.py`` from your ``PYTHONPATH``; see
1926  ``baseDebug`` for more about ``debug.py`` files.
1927 
1928  This task supports the following debug variables:
1929 
1930  - ``saveCountIm``
1931  If True then save the Epoch Count Image as a fits file in the `figPath`
1932  - ``figPath``
1933  Path to save the debug fits images and figures
1934 
1935  For example, put something like:
1936 
1937  .. code-block:: python
1938 
1939  import lsstDebug
1940  def DebugInfo(name):
1941  di = lsstDebug.getInfo(name)
1942  if name == "lsst.pipe.tasks.assembleCoadd":
1943  di.saveCountIm = True
1944  di.figPath = "/desired/path/to/debugging/output/images"
1945  return di
1946  lsstDebug.Info = DebugInfo
1947 
1948  into your ``debug.py`` file and run ``assemebleCoadd.py`` with the
1949  ``--debug`` flag. Some subtasks may have their own debug variables;
1950  see individual Task documentation.
1951 
1952  Examples
1953  --------
1954  ``CompareWarpAssembleCoaddTask`` assembles a set of warped images into a
1955  coadded image. The ``CompareWarpAssembleCoaddTask`` is invoked by running
1956  ``assembleCoadd.py`` with the flag ``--compareWarpCoadd``.
1957  Usage of ``assembleCoadd.py`` expects a data reference to the tract patch
1958  and filter to be coadded (specified using
1959  '--id = [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]')
1960  along with a list of coaddTempExps to attempt to coadd (specified using
1961  '--selectId [KEY=VALUE1[^VALUE2[^VALUE3...] [KEY=VALUE1[^VALUE2[^VALUE3...] ...]]').
1962  Only the warps that cover the specified tract and patch will be coadded.
1963  A list of the available optional arguments can be obtained by calling
1964  ``assembleCoadd.py`` with the ``--help`` command line argument:
1965 
1966  .. code-block:: none
1967 
1968  assembleCoadd.py --help
1969 
1970  To demonstrate usage of the ``CompareWarpAssembleCoaddTask`` in the larger
1971  context of multi-band processing, we will generate the HSC-I & -R band
1972  oadds from HSC engineering test data provided in the ``ci_hsc`` package.
1973  To begin, assuming that the lsst stack has been already set up, we must
1974  set up the ``obs_subaru`` and ``ci_hsc`` packages.
1975  This defines the environment variable ``$CI_HSC_DIR`` and points at the
1976  location of the package. The raw HSC data live in the ``$CI_HSC_DIR/raw``
1977  directory. To begin assembling the coadds, we must first
1978 
1979  - processCcd
1980  process the individual ccds in $CI_HSC_RAW to produce calibrated exposures
1981  - makeSkyMap
1982  create a skymap that covers the area of the sky present in the raw exposures
1983  - makeCoaddTempExp
1984  warp the individual calibrated exposures to the tangent plane of the coadd
1985 
1986  We can perform all of these steps by running
1987 
1988  .. code-block:: none
1989 
1990  $CI_HSC_DIR scons warp-903986 warp-904014 warp-903990 warp-904010 warp-903988
1991 
1992  This will produce warped ``coaddTempExps`` for each visit. To coadd the
1993  warped data, we call ``assembleCoadd.py`` as follows:
1994 
1995  .. code-block:: none
1996 
1997  assembleCoadd.py --compareWarpCoadd $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-I \
1998  --selectId visit=903986 ccd=16 --selectId visit=903986 ccd=22 --selectId visit=903986 ccd=23 \
1999  --selectId visit=903986 ccd=100 --selectId visit=904014 ccd=1 --selectId visit=904014 ccd=6 \
2000  --selectId visit=904014 ccd=12 --selectId visit=903990 ccd=18 --selectId visit=903990 ccd=25 \
2001  --selectId visit=904010 ccd=4 --selectId visit=904010 ccd=10 --selectId visit=904010 ccd=100 \
2002  --selectId visit=903988 ccd=16 --selectId visit=903988 ccd=17 --selectId visit=903988 ccd=23 \
2003  --selectId visit=903988 ccd=24
2004 
2005  This will process the HSC-I band data. The results are written in
2006  ``$CI_HSC_DIR/DATA/deepCoadd-results/HSC-I``.
2007  """
2008  ConfigClass = CompareWarpAssembleCoaddConfig
2009  _DefaultName = "compareWarpAssembleCoadd"
2010 
2011  def __init__(self, *args, **kwargs):
2012  AssembleCoaddTask.__init__(self, *args, **kwargs)
2013  self.makeSubtask("assembleStaticSkyModel")
2014  detectionSchema = afwTable.SourceTable.makeMinimalSchema()
2015  self.makeSubtask("detect", schema=detectionSchema)
2016  if self.config.doPreserveContainedBySource:
2017  self.makeSubtask("detectTemplate", schema=afwTable.SourceTable.makeMinimalSchema())
2018  if self.config.doScaleWarpVariance:
2019  self.makeSubtask("scaleWarpVariance")
2020  if self.config.doFilterMorphological:
2021  self.makeSubtask("maskStreaks")
2022 
2023  @utils.inheritDoc(AssembleCoaddTask)
2024  def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs):
2025  """
2026  Generate a templateCoadd to use as a naive model of static sky to
2027  subtract from PSF-Matched warps.
2028 
2029  Returns
2030  -------
2031  result : `lsst.pipe.base.Struct`
2032  Result struct with components:
2033 
2034  - ``templateCoadd`` : coadded exposure (``lsst.afw.image.Exposure``)
2035  - ``nImage`` : N Image (``lsst.afw.image.Image``)
2036  """
2037  # Ensure that psfMatchedWarps are used as input warps for template generation
2038  staticSkyModelInputRefs = copy.deepcopy(inputRefs)
2039  staticSkyModelInputRefs.inputWarps = inputRefs.psfMatchedWarps
2040 
2041  # Because subtasks don't have connections we have to make one.
2042  # The main task's `templateCoadd` is the subtask's `coaddExposure`
2043  staticSkyModelOutputRefs = copy.deepcopy(outputRefs)
2044  if self.config.assembleStaticSkyModel.doWrite:
2045  staticSkyModelOutputRefs.coaddExposure = staticSkyModelOutputRefs.templateCoadd
2046  # Remove template coadd from both subtask's and main tasks outputs,
2047  # because it is handled by the subtask as `coaddExposure`
2048  del outputRefs.templateCoadd
2049  del staticSkyModelOutputRefs.templateCoadd
2050 
2051  # A PSF-Matched nImage does not exist as a dataset type
2052  if 'nImage' in staticSkyModelOutputRefs.keys():
2053  del staticSkyModelOutputRefs.nImage
2054 
2055  templateCoadd = self.assembleStaticSkyModel.runQuantum(butlerQC, staticSkyModelInputRefs,
2056  staticSkyModelOutputRefs)
2057  if templateCoadd is None:
2058  raise RuntimeError(self._noTemplateMessage(self.assembleStaticSkyModel.warpType))
2059 
2060  return pipeBase.Struct(templateCoadd=templateCoadd.coaddExposure,
2061  nImage=templateCoadd.nImage,
2062  warpRefList=templateCoadd.warpRefList,
2063  imageScalerList=templateCoadd.imageScalerList,
2064  weightList=templateCoadd.weightList)
2065 
2066  @utils.inheritDoc(AssembleCoaddTask)
2067  def makeSupplementaryData(self, dataRef, selectDataList=None, warpRefList=None):
2068  """
2069  Generate a templateCoadd to use as a naive model of static sky to
2070  subtract from PSF-Matched warps.
2071 
2072  Returns
2073  -------
2074  result : `lsst.pipe.base.Struct`
2075  Result struct with components:
2076 
2077  - ``templateCoadd``: coadded exposure (``lsst.afw.image.Exposure``)
2078  - ``nImage``: N Image (``lsst.afw.image.Image``)
2079  """
2080  templateCoadd = self.assembleStaticSkyModel.runDataRef(dataRef, selectDataList, warpRefList)
2081  if templateCoadd is None:
2082  raise RuntimeError(self._noTemplateMessage(self.assembleStaticSkyModel.warpType))
2083 
2084  return pipeBase.Struct(templateCoadd=templateCoadd.coaddExposure,
2085  nImage=templateCoadd.nImage,
2086  warpRefList=templateCoadd.warpRefList,
2087  imageScalerList=templateCoadd.imageScalerList,
2088  weightList=templateCoadd.weightList)
2089 
2090  def _noTemplateMessage(self, warpType):
2091  warpName = (warpType[0].upper() + warpType[1:])
2092  message = """No %(warpName)s warps were found to build the template coadd which is
2093  required to run CompareWarpAssembleCoaddTask. To continue assembling this type of coadd,
2094  first either rerun makeCoaddTempExp with config.make%(warpName)s=True or
2095  coaddDriver with config.makeCoadTempExp.make%(warpName)s=True, before assembleCoadd.
2096 
2097  Alternatively, to use another algorithm with existing warps, retarget the CoaddDriverConfig to
2098  another algorithm like:
2099 
2100  from lsst.pipe.tasks.assembleCoadd import SafeClipAssembleCoaddTask
2101  config.assemble.retarget(SafeClipAssembleCoaddTask)
2102  """ % {"warpName": warpName}
2103  return message
2104 
2105  @utils.inheritDoc(AssembleCoaddTask)
2106  @pipeBase.timeMethod
2107  def run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
2108  supplementaryData, *args, **kwargs):
2109  """Assemble the coadd.
2110 
2111  Find artifacts and apply them to the warps' masks creating a list of
2112  alternative masks with a new "CLIPPED" plane and updated "NO_DATA"
2113  plane. Then pass these alternative masks to the base class's `run`
2114  method.
2115 
2116  The input parameters ``supplementaryData`` is a `lsst.pipe.base.Struct`
2117  that must contain a ``templateCoadd`` that serves as the
2118  model of the static sky.
2119  """
2120 
2121  # Check and match the order of the supplementaryData
2122  # (PSF-matched) inputs to the order of the direct inputs,
2123  # so that the artifact mask is applied to the right warp
2124  dataIds = [ref.dataId for ref in tempExpRefList]
2125  psfMatchedDataIds = [ref.dataId for ref in supplementaryData.warpRefList]
2126 
2127  if dataIds != psfMatchedDataIds:
2128  self.log.info("Reordering and or/padding PSF-matched visit input list")
2129  supplementaryData.warpRefList = reorderAndPadList(supplementaryData.warpRefList,
2130  psfMatchedDataIds, dataIds)
2131  supplementaryData.imageScalerList = reorderAndPadList(supplementaryData.imageScalerList,
2132  psfMatchedDataIds, dataIds)
2133 
2134  # Use PSF-Matched Warps (and corresponding scalers) and coadd to find artifacts
2135  spanSetMaskList = self.findArtifacts(supplementaryData.templateCoadd,
2136  supplementaryData.warpRefList,
2137  supplementaryData.imageScalerList)
2138 
2139  badMaskPlanes = self.config.badMaskPlanes[:]
2140  badMaskPlanes.append("CLIPPED")
2141  badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes)
2142 
2143  result = AssembleCoaddTask.run(self, skyInfo, tempExpRefList, imageScalerList, weightList,
2144  spanSetMaskList, mask=badPixelMask)
2145 
2146  # Propagate PSF-matched EDGE pixels to coadd SENSOR_EDGE and INEXACT_PSF
2147  # Psf-Matching moves the real edge inwards
2148  self.applyAltEdgeMask(result.coaddExposure.maskedImage.mask, spanSetMaskList)
2149  return result
2150 
2151  def applyAltEdgeMask(self, mask, altMaskList):
2152  """Propagate alt EDGE mask to SENSOR_EDGE AND INEXACT_PSF planes.
2153 
2154  Parameters
2155  ----------
2156  mask : `lsst.afw.image.Mask`
2157  Original mask.
2158  altMaskList : `list`
2159  List of Dicts containing ``spanSet`` lists.
2160  Each element contains the new mask plane name (e.g. "CLIPPED
2161  and/or "NO_DATA") as the key, and list of ``SpanSets`` to apply to
2162  the mask.
2163  """
2164  maskValue = mask.getPlaneBitMask(["SENSOR_EDGE", "INEXACT_PSF"])
2165  for visitMask in altMaskList:
2166  if "EDGE" in visitMask:
2167  for spanSet in visitMask['EDGE']:
2168  spanSet.clippedTo(mask.getBBox()).setMask(mask, maskValue)
2169 
2170  def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList):
2171  """Find artifacts.
2172 
2173  Loop through warps twice. The first loop builds a map with the count
2174  of how many epochs each pixel deviates from the templateCoadd by more
2175  than ``config.chiThreshold`` sigma. The second loop takes each
2176  difference image and filters the artifacts detected in each using
2177  count map to filter out variable sources and sources that are
2178  difficult to subtract cleanly.
2179 
2180  Parameters
2181  ----------
2182  templateCoadd : `lsst.afw.image.Exposure`
2183  Exposure to serve as model of static sky.
2184  tempExpRefList : `list`
2185  List of data references to warps.
2186  imageScalerList : `list`
2187  List of image scalers.
2188 
2189  Returns
2190  -------
2191  altMasks : `list`
2192  List of dicts containing information about CLIPPED
2193  (i.e., artifacts), NO_DATA, and EDGE pixels.
2194  """
2195 
2196  self.log.debug("Generating Count Image, and mask lists.")
2197  coaddBBox = templateCoadd.getBBox()
2198  slateIm = afwImage.ImageU(coaddBBox)
2199  epochCountImage = afwImage.ImageU(coaddBBox)
2200  nImage = afwImage.ImageU(coaddBBox)
2201  spanSetArtifactList = []
2202  spanSetNoDataMaskList = []
2203  spanSetEdgeList = []
2204  spanSetBadMorphoList = []
2205  badPixelMask = self.getBadPixelMask()
2206 
2207  # mask of the warp diffs should = that of only the warp
2208  templateCoadd.mask.clearAllMaskPlanes()
2209 
2210  if self.config.doPreserveContainedBySource:
2211  templateFootprints = self.detectTemplate.detectFootprints(templateCoadd)
2212  else:
2213  templateFootprints = None
2214 
2215  for warpRef, imageScaler in zip(tempExpRefList, imageScalerList):
2216  warpDiffExp = self._readAndComputeWarpDiff(warpRef, imageScaler, templateCoadd)
2217  if warpDiffExp is not None:
2218  # This nImage only approximates the final nImage because it uses the PSF-matched mask
2219  nImage.array += (numpy.isfinite(warpDiffExp.image.array)
2220  * ((warpDiffExp.mask.array & badPixelMask) == 0)).astype(numpy.uint16)
2221  fpSet = self.detect.detectFootprints(warpDiffExp, doSmooth=False, clearMask=True)
2222  fpSet.positive.merge(fpSet.negative)
2223  footprints = fpSet.positive
2224  slateIm.set(0)
2225  spanSetList = [footprint.spans for footprint in footprints.getFootprints()]
2226 
2227  # Remove artifacts due to defects before they contribute to the epochCountImage
2228  if self.config.doPrefilterArtifacts:
2229  spanSetList = self.prefilterArtifacts(spanSetList, warpDiffExp)
2230 
2231  # Clear mask before adding prefiltered spanSets
2232  self.detect.clearMask(warpDiffExp.mask)
2233  for spans in spanSetList:
2234  spans.setImage(slateIm, 1, doClip=True)
2235  spans.setMask(warpDiffExp.mask, warpDiffExp.mask.getPlaneBitMask("DETECTED"))
2236  epochCountImage += slateIm
2237 
2238  if self.config.doFilterMorphological:
2239  maskName = self.config.streakMaskName
2240  _ = self.maskStreaks.run(warpDiffExp)
2241  streakMask = warpDiffExp.mask
2242  spanSetStreak = afwGeom.SpanSet.fromMask(streakMask,
2243  streakMask.getPlaneBitMask(maskName)).split()
2244 
2245  # PSF-Matched warps have less available area (~the matching kernel) because the calexps
2246  # undergo a second convolution. Pixels with data in the direct warp
2247  # but not in the PSF-matched warp will not have their artifacts detected.
2248  # NaNs from the PSF-matched warp therefore must be masked in the direct warp
2249  nans = numpy.where(numpy.isnan(warpDiffExp.maskedImage.image.array), 1, 0)
2250  nansMask = afwImage.makeMaskFromArray(nans.astype(afwImage.MaskPixel))
2251  nansMask.setXY0(warpDiffExp.getXY0())
2252  edgeMask = warpDiffExp.mask
2253  spanSetEdgeMask = afwGeom.SpanSet.fromMask(edgeMask,
2254  edgeMask.getPlaneBitMask("EDGE")).split()
2255  else:
2256  # If the directWarp has <1% coverage, the psfMatchedWarp can have 0% and not exist
2257  # In this case, mask the whole epoch
2258  nansMask = afwImage.MaskX(coaddBBox, 1)
2259  spanSetList = []
2260  spanSetEdgeMask = []
2261  spanSetStreak = []
2262 
2263  spanSetNoDataMask = afwGeom.SpanSet.fromMask(nansMask).split()
2264 
2265  spanSetNoDataMaskList.append(spanSetNoDataMask)
2266  spanSetArtifactList.append(spanSetList)
2267  spanSetEdgeList.append(spanSetEdgeMask)
2268  if self.config.doFilterMorphological:
2269  spanSetBadMorphoList.append(spanSetStreak)
2270 
2271  if lsstDebug.Info(__name__).saveCountIm:
2272  path = self._dataRef2DebugPath("epochCountIm", tempExpRefList[0], coaddLevel=True)
2273  epochCountImage.writeFits(path)
2274 
2275  for i, spanSetList in enumerate(spanSetArtifactList):
2276  if spanSetList:
2277  filteredSpanSetList = self.filterArtifacts(spanSetList, epochCountImage, nImage,
2278  templateFootprints)
2279  spanSetArtifactList[i] = filteredSpanSetList
2280  if self.config.doFilterMorphological:
2281  spanSetArtifactList[i] += spanSetBadMorphoList[i]
2282 
2283  altMasks = []
2284  for artifacts, noData, edge in zip(spanSetArtifactList, spanSetNoDataMaskList, spanSetEdgeList):
2285  altMasks.append({'CLIPPED': artifacts,
2286  'NO_DATA': noData,
2287  'EDGE': edge})
2288  return altMasks
2289 
2290  def prefilterArtifacts(self, spanSetList, exp):
2291  """Remove artifact candidates covered by bad mask plane.
2292 
2293  Any future editing of the candidate list that does not depend on
2294  temporal information should go in this method.
2295 
2296  Parameters
2297  ----------
2298  spanSetList : `list`
2299  List of SpanSets representing artifact candidates.
2300  exp : `lsst.afw.image.Exposure`
2301  Exposure containing mask planes used to prefilter.
2302 
2303  Returns
2304  -------
2305  returnSpanSetList : `list`
2306  List of SpanSets with artifacts.
2307  """
2308  badPixelMask = exp.mask.getPlaneBitMask(self.config.prefilterArtifactsMaskPlanes)
2309  goodArr = (exp.mask.array & badPixelMask) == 0
2310  returnSpanSetList = []
2311  bbox = exp.getBBox()
2312  x0, y0 = exp.getXY0()
2313  for i, span in enumerate(spanSetList):
2314  y, x = span.clippedTo(bbox).indices()
2315  yIndexLocal = numpy.array(y) - y0
2316  xIndexLocal = numpy.array(x) - x0
2317  goodRatio = numpy.count_nonzero(goodArr[yIndexLocal, xIndexLocal])/span.getArea()
2318  if goodRatio > self.config.prefilterArtifactsRatio:
2319  returnSpanSetList.append(span)
2320  return returnSpanSetList
2321 
2322  def filterArtifacts(self, spanSetList, epochCountImage, nImage, footprintsToExclude=None):
2323  """Filter artifact candidates.
2324 
2325  Parameters
2326  ----------
2327  spanSetList : `list`
2328  List of SpanSets representing artifact candidates.
2329  epochCountImage : `lsst.afw.image.Image`
2330  Image of accumulated number of warpDiff detections.
2331  nImage : `lsst.afw.image.Image`
2332  Image of the accumulated number of total epochs contributing.
2333 
2334  Returns
2335  -------
2336  maskSpanSetList : `list`
2337  List of SpanSets with artifacts.
2338  """
2339 
2340  maskSpanSetList = []
2341  x0, y0 = epochCountImage.getXY0()
2342  for i, span in enumerate(spanSetList):
2343  y, x = span.indices()
2344  yIdxLocal = [y1 - y0 for y1 in y]
2345  xIdxLocal = [x1 - x0 for x1 in x]
2346  outlierN = epochCountImage.array[yIdxLocal, xIdxLocal]
2347  totalN = nImage.array[yIdxLocal, xIdxLocal]
2348 
2349  # effectiveMaxNumEpochs is broken line (fraction of N) with characteristic config.maxNumEpochs
2350  effMaxNumEpochsHighN = (self.config.maxNumEpochs
2351  + self.config.maxFractionEpochsHigh*numpy.mean(totalN))
2352  effMaxNumEpochsLowN = self.config.maxFractionEpochsLow * numpy.mean(totalN)
2353  effectiveMaxNumEpochs = int(min(effMaxNumEpochsLowN, effMaxNumEpochsHighN))
2354  nPixelsBelowThreshold = numpy.count_nonzero((outlierN > 0)
2355  & (outlierN <= effectiveMaxNumEpochs))
2356  percentBelowThreshold = nPixelsBelowThreshold / len(outlierN)
2357  if percentBelowThreshold > self.config.spatialThreshold:
2358  maskSpanSetList.append(span)
2359 
2360  if self.config.doPreserveContainedBySource and footprintsToExclude is not None:
2361  # If a candidate is contained by a footprint on the template coadd, do not clip
2362  filteredMaskSpanSetList = []
2363  for span in maskSpanSetList:
2364  doKeep = True
2365  for footprint in footprintsToExclude.positive.getFootprints():
2366  if footprint.spans.contains(span):
2367  doKeep = False
2368  break
2369  if doKeep:
2370  filteredMaskSpanSetList.append(span)
2371  maskSpanSetList = filteredMaskSpanSetList
2372 
2373  return maskSpanSetList
2374 
2375  def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd):
2376  """Fetch a warp from the butler and return a warpDiff.
2377 
2378  Parameters
2379  ----------
2380  warpRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
2381  Butler dataRef for the warp.
2382  imageScaler : `lsst.pipe.tasks.scaleZeroPoint.ImageScaler`
2383  An image scaler object.
2384  templateCoadd : `lsst.afw.image.Exposure`
2385  Exposure to be substracted from the scaled warp.
2386 
2387  Returns
2388  -------
2389  warp : `lsst.afw.image.Exposure`
2390  Exposure of the image difference between the warp and template.
2391  """
2392 
2393  # If the PSF-Matched warp did not exist for this direct warp
2394  # None is holding its place to maintain order in Gen 3
2395  if warpRef is None:
2396  return None
2397  # Warp comparison must use PSF-Matched Warps regardless of requested coadd warp type
2398  warpName = self.getTempExpDatasetName('psfMatched')
2399  if not isinstance(warpRef, DeferredDatasetHandle):
2400  if not warpRef.datasetExists(warpName):
2401  self.log.warn("Could not find %s %s; skipping it", warpName, warpRef.dataId)
2402  return None
2403  warp = warpRef.get(datasetType=warpName, immediate=True)
2404  # direct image scaler OK for PSF-matched Warp
2405  imageScaler.scaleMaskedImage(warp.getMaskedImage())
2406  mi = warp.getMaskedImage()
2407  if self.config.doScaleWarpVariance:
2408  try:
2409  self.scaleWarpVariance.run(mi)
2410  except Exception as exc:
2411  self.log.warn("Unable to rescale variance of warp (%s); leaving it as-is" % (exc,))
2412  mi -= templateCoadd.getMaskedImage()
2413  return warp
2414 
2415  def _dataRef2DebugPath(self, prefix, warpRef, coaddLevel=False):
2416  """Return a path to which to write debugging output.
2417 
2418  Creates a hyphen-delimited string of dataId values for simple filenames.
2419 
2420  Parameters
2421  ----------
2422  prefix : `str`
2423  Prefix for filename.
2424  warpRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
2425  Butler dataRef to make the path from.
2426  coaddLevel : `bool`, optional.
2427  If True, include only coadd-level keys (e.g., 'tract', 'patch',
2428  'filter', but no 'visit').
2429 
2430  Returns
2431  -------
2432  result : `str`
2433  Path for debugging output.
2434  """
2435  if coaddLevel:
2436  keys = warpRef.getButler().getKeys(self.getCoaddDatasetName(self.warpType))
2437  else:
2438  keys = warpRef.dataId.keys()
2439  keyList = sorted(keys, reverse=True)
2440  directory = lsstDebug.Info(__name__).figPath if lsstDebug.Info(__name__).figPath else "."
2441  filename = "%s-%s.fits" % (prefix, '-'.join([str(warpRef.dataId[k]) for k in keyList]))
2442  return os.path.join(directory, filename)
2443 
2444 
2445 def reorderAndPadList(inputList, inputKeys, outputKeys, padWith=None):
2446  """Match the order of one list to another, padding if necessary
2447 
2448  Parameters
2449  ----------
2450  inputList : list
2451  List to be reordered and padded. Elements can be any type.
2452  inputKeys : iterable
2453  Iterable of values to be compared with outputKeys.
2454  Length must match `inputList`
2455  outputKeys : iterable
2456  Iterable of values to be compared with inputKeys.
2457  padWith :
2458  Any value to be inserted where inputKey not in outputKeys
2459 
2460  Returns
2461  -------
2462  list
2463  Copy of inputList reordered per outputKeys and padded with `padWith`
2464  so that the length matches length of outputKeys.
2465  """
2466  outputList = []
2467  for d in outputKeys:
2468  if d in inputKeys:
2469  outputList.append(inputList[inputKeys.index(d)])
2470  else:
2471  outputList.append(padWith)
2472  return outputList
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask.filterArtifacts
def filterArtifacts(self, spanSetList, epochCountImage, nImage, footprintsToExclude=None)
Definition: assembleCoadd.py:2322
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::log.log.logContinued.warn
def warn(fmt, *args)
Definition: logContinued.py:205
lsst.pipe.tasks.assembleCoadd.getTempExpRefList
def getTempExpRefList(self, patchRef, calExpRefList)
Definition: assembleCoadd.py:588
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig
Definition: assembleCoadd.py:1740
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig
Definition: assembleCoadd.py:1231
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask.makeSupplementaryData
def makeSupplementaryData(self, dataRef, selectDataList=None, warpRefList=None)
Definition: assembleCoadd.py:2067
lsst::log.log.logContinued.info
def info(fmt, *args)
Definition: logContinued.py:201
lsst::afw::image::Mask
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
lsst.pex.config.wrap.setDefaults
setDefaults
Definition: wrap.py:293
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddTask.run
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, *args, **kwargs)
Definition: assembleCoadd.py:1434
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask.run
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, supplementaryData, *args, **kwargs)
Definition: assembleCoadd.py:2107
lsst.pipe.tasks.assembleCoadd.assembleSubregion
def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList, weightList, altMaskList, statsFlags, statsCtrl, nImage=None)
Definition: assembleCoadd.py:862
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig.validate
def validate(self)
Definition: assembleCoadd.py:1874
lsst.gdb.afw.printers.debug
bool debug
Definition: printers.py:9
ast::append
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddTask.__init__
def __init__(self, *args, **kwargs)
Definition: assembleCoadd.py:1425
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddTask
Definition: assembleCoadd.py:1313
lsst.pipe.tasks.assembleCoadd.processResults
def processResults(self, coaddExposure, brightObjectMasks=None, dataId=None)
Definition: assembleCoadd.py:525
lsst::afw::math::makeStatistics
Statistics makeStatistics(lsst::afw::math::MaskedVector< EntryT > const &mv, std::vector< WeightPixel > const &vweights, int const flags, StatisticsControl const &sctrl=StatisticsControl())
The makeStatistics() overload to handle lsst::afw::math::MaskedVector<>
Definition: Statistics.h:520
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.clipDetection
clipDetection
Definition: assembleCoadd.py:1242
astshim.keyMap.keyMapContinued.keys
def keys(self)
Definition: keyMapContinued.py:6
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask.findArtifacts
def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList)
Definition: assembleCoadd.py:2170
lsst::afw::math::stringToStatisticsProperty
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
Definition: Statistics.cc:747
lsst.pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
lsst.pipe.tasks.coaddBase.makeSkyInfo
def makeSkyInfo(skyMap, tractId, patchId)
Definition: coaddBase.py:279
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig.statistic
statistic
Definition: assembleCoadd.py:1842
lsst.pipe.tasks.assembleCoadd.AssembleCoaddConnections
Definition: assembleCoadd.py:56
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask._noTemplateMessage
def _noTemplateMessage(self, warpType)
Definition: assembleCoadd.py:2090
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig.detectTemplate
detectTemplate
Definition: assembleCoadd.py:1750
lsst.pipe.tasks.assembleCoadd.reorderAndPadList
def reorderAndPadList(inputList, inputKeys, outputKeys, padWith=None)
Definition: assembleCoadd.py:2445
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig.warpType
warpType
Definition: assembleCoadd.py:1880
lsst::geom::all
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
Definition: CoordinateExpr.h:81
lsst.pipe.tasks.assembleCoadd.removeMaskPlanes
def removeMaskPlanes(self, maskedImage)
Definition: assembleCoadd.py:939
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddTask.detectClip
def detectClip(self, exp, tempExpRefList)
Definition: assembleCoadd.py:1526
lsst.pipe.tasks.assembleCoadd.assembleMetadata
def assembleMetadata(self, coaddExposure, tempExpRefList, weightList)
Definition: assembleCoadd.py:801
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.sigmaClip
sigmaClip
Definition: assembleCoadd.py:1293
lsst::afw::detection::FootprintSet
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.assembleMeanClipCoadd
assembleMeanClipCoadd
Definition: assembleCoadd.py:1238
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask
Definition: assembleCoadd.py:1887
lsst::afw::geom::SpanSet
A compact representation of a collection of pixels.
Definition: SpanSet.h:78
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConnections.__init__
def __init__(self, *config=None)
Definition: assembleCoadd.py:1732
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask.prefilterArtifacts
def prefilterArtifacts(self, spanSetList, exp)
Definition: assembleCoadd.py:2290
lsst::log.log.logContinued.fatal
def fatal(fmt, *args)
Definition: logContinued.py:217
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask._readAndComputeWarpDiff
def _readAndComputeWarpDiff(self, warpRef, imageScaler, templateCoadd)
Definition: assembleCoadd.py:2375
lsst.pipe.tasks.assembleCoadd.run
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)
Definition: assembleCoadd.py:720
lsst.pipe.tasks.assembleCoadd.prepareStats
def prepareStats(self, mask=None)
Definition: assembleCoadd.py:687
lsstDebug.Info
Definition: lsstDebug.py:28
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig.detect
detect
Definition: assembleCoadd.py:1746
lsst.pex.config
Definition: __init__.py:1
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig.doUsePsfMatchedPolygons
doUsePsfMatchedPolygons
Definition: assembleCoadd.py:1843
lsst.pipe.tasks.assembleCoadd.makeSupplementaryData
def makeSupplementaryData(self, dataRef, selectDataList=None, warpRefList=None)
Definition: assembleCoadd.py:545
lsst::afw::detection::Threshold
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:43
lsst::log
Definition: Log.h:706
max
int max
Definition: BoundedField.cc:104
lsst::afw::table
Definition: table.dox:3
std::int32_t
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConnections
Definition: assembleCoadd.py:1713
lsst::utils
Definition: Backtrace.h:29
lsst.pipe.tasks.assembleCoadd.setInexactPsf
def setInexactPsf(self, mask)
Definition: assembleCoadd.py:1102
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.assembleMeanCoadd
assembleMeanCoadd
Definition: assembleCoadd.py:1234
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddTask.buildDifferenceImage
def buildDifferenceImage(self, skyInfo, tempExpRefList, imageScalerList, weightList)
Definition: assembleCoadd.py:1490
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.doSigmaClip
doSigmaClip
Definition: assembleCoadd.py:1305
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddTask.detectClipBig
def detectClipBig(self, clipList, clipFootprints, clipIndices, detectionFootprints, maskClipValue, maskDetValue, coaddBBox)
Definition: assembleCoadd.py:1652
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask.applyAltEdgeMask
def applyAltEdgeMask(self, mask, altMaskList)
Definition: assembleCoadd.py:2151
lsst::afw::detection
Definition: Footprint.h:50
lsst.pipe.tasks.coaddHelpers.getGroupDataRef
def getGroupDataRef(butler, datasetType, groupTuple, keys)
Definition: coaddHelpers.py:99
lsst::afw::math::StatisticsControl
Pass parameters to a Statistics object.
Definition: Statistics.h:93
lsst::geom
Definition: AffineTransform.h:36
lsst.pipe.tasks.assembleCoadd.countMaskFromFootprint
def countMaskFromFootprint(mask, footprint, bitmask, ignoreMask)
Definition: assembleCoadd.py:1198
lsst.pipe.tasks.assembleCoadd.readBrightObjectMasks
def readBrightObjectMasks(self, dataRef)
Definition: assembleCoadd.py:1038
items
std::vector< SchemaItem< Flag > > * items
Definition: BaseColumnView.cc:142
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.setDefaults
def setDefaults(self)
Definition: assembleCoadd.py:1275
lsst.pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
lsst.pipe.tasks.assembleCoadd.AssembleCoaddDataIdContainer
Definition: assembleCoadd.py:1170
lsst::afw::image.slicing.Factory
Factory
Definition: slicing.py:252
list
daf::base::PropertyList * list
Definition: fits.cc:913
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.validate
def validate(self)
Definition: assembleCoadd.py:1301
min
int min
Definition: BoundedField.cc:103
lsst.pipe.tasks.assembleCoadd.prepareInputs
def prepareInputs(self, refList)
Definition: assembleCoadd.py:612
lsst::afw::math
Definition: statistics.dox:6
lsst.pipe.tasks.coaddBase.makeCoaddSuffix
def makeCoaddSuffix(warpType="direct")
Definition: coaddBase.py:336
lsst.pipe.tasks.assembleCoadd.makeSupplementaryDataGen3
def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs)
Definition: assembleCoadd.py:564
lsst::geom::Point< int, 2 >
lsst.pipe.tasks.coaddBase.CoaddBaseTask
Base class for coaddition.
Definition: coaddBase.py:131
lsst.pipe.tasks.assembleCoadd.setBrightObjectMasks
def setBrightObjectMasks(self, exposure, brightObjectMasks, dataId=None)
Definition: assembleCoadd.py:1060
lsst::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:55
lsst::coadd::utils::setCoaddEdgeBits
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
Definition: setCoaddEdgeBits.cc:42
lsst.pex::exceptions
Definition: Exception.h:37
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask._dataRef2DebugPath
def _dataRef2DebugPath(self, prefix, warpRef, coaddLevel=False)
Definition: assembleCoadd.py:2415
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.clipIter
clipIter
Definition: assembleCoadd.py:1294
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask.__init__
def __init__(self, *args, **kwargs)
Definition: assembleCoadd.py:2011
lsst.pipe.tasks.assembleCoadd.SafeClipAssembleCoaddConfig.statistic
statistic
Definition: assembleCoadd.py:1295
lsst::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
lsst.pipe.base
Definition: __init__.py:1
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig.assembleStaticSkyModel
assembleStaticSkyModel
Definition: assembleCoadd.py:1741
lsst::afw::math::statisticsStack
std::shared_ptr< image::MaskedImage< PixelT > > statisticsStack(image::MaskedImage< PixelT > const &image, Property flags, char dimension, StatisticsControl const &sctrl)
A function to compute statistics on the rows or columns of an image.
Definition: Stack.cc:488
lsst.pipe.tasks.assembleCoadd.applyAltMaskPlanes
def applyAltMaskPlanes(self, mask, altMaskSpans)
Definition: assembleCoadd.py:984
lsst::meas::algorithms
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
Definition: CoaddBoundedField.h:34
lsst.pex.config.wrap.validate
validate
Definition: wrap.py:295
lsst.pipe.tasks.coaddHelpers.groupPatchExposures
def groupPatchExposures(patchDataRef, calexpDataRefList, coaddDatasetType="deepCoadd", tempExpDatasetType="deepCoadd_directWarp")
Definition: coaddHelpers.py:59
lsst::geom::Extent< int, 2 >
lsst::coadd::utils
Definition: addToCoadd.h:37
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddConfig.setDefaults
def setDefaults(self)
Definition: assembleCoadd.py:1840
lsst.pipe.tasks.assembleCoadd.AssembleCoaddDataIdContainer.makeDataRefList
def makeDataRefList(self, namespace)
Definition: assembleCoadd.py:1174
lsst.pipe.tasks.assembleCoadd.shrinkValidPolygons
def shrinkValidPolygons(self, coaddInputs)
Definition: assembleCoadd.py:1016
lsst.pipe.tasks.assembleCoadd.CompareWarpAssembleCoaddTask.makeSupplementaryDataGen3
def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs)
Definition: assembleCoadd.py:2024
lsst::afw::geom
Definition: frameSetUtils.h:40