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