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