LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
dcrAssembleCoadd.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 
23 from math import ceil
24 import numpy as np
25 from scipy import ndimage
26 import lsst.geom as geom
27 import lsst.afw.image as afwImage
28 import lsst.afw.table as afwTable
29 import lsst.coadd.utils as coaddUtils
30 from lsst.daf.butler import DeferredDatasetHandle
31 from lsst.ip.diffim.dcrModel import applyDcr, calculateDcr, DcrModel
32 import lsst.meas.algorithms as measAlg
33 from lsst.meas.base import SingleFrameMeasurementTask
34 import lsst.pex.config as pexConfig
35 import lsst.pipe.base as pipeBase
36 import lsst.utils as utils
37 from lsst.utils.timer import timeMethod
38 from .assembleCoadd import (AssembleCoaddConnections,
39  AssembleCoaddTask,
40  CompareWarpAssembleCoaddConfig,
41  CompareWarpAssembleCoaddTask)
42 from .coaddBase import makeSkyInfo
43 from .measurePsf import MeasurePsfTask
44 
45 __all__ = ["DcrAssembleCoaddConnections", "DcrAssembleCoaddTask", "DcrAssembleCoaddConfig"]
46 
47 
49  dimensions=("tract", "patch", "band", "skymap"),
50  defaultTemplates={"inputWarpName": "deep",
51  "inputCoaddName": "deep",
52  "outputCoaddName": "dcr",
53  "warpType": "direct",
54  "warpTypeSuffix": "",
55  "fakesType": ""}):
56  inputWarps = pipeBase.connectionTypes.Input(
57  doc=("Input list of warps to be assembled i.e. stacked."
58  "Note that this will often be different than the inputCoaddName."
59  "WarpType (e.g. direct, psfMatched) is controlled by the warpType config parameter"),
60  name="{inputWarpName}Coadd_{warpType}Warp",
61  storageClass="ExposureF",
62  dimensions=("tract", "patch", "skymap", "visit", "instrument"),
63  deferLoad=True,
64  multiple=True
65  )
66  templateExposure = pipeBase.connectionTypes.Input(
67  doc="Input coadded exposure, produced by previous call to AssembleCoadd",
68  name="{fakesType}{inputCoaddName}Coadd{warpTypeSuffix}",
69  storageClass="ExposureF",
70  dimensions=("tract", "patch", "skymap", "band"),
71  )
72  dcrCoadds = pipeBase.connectionTypes.Output(
73  doc="Output coadded exposure, produced by stacking input warps",
74  name="{fakesType}{outputCoaddName}Coadd{warpTypeSuffix}",
75  storageClass="ExposureF",
76  dimensions=("tract", "patch", "skymap", "band", "subfilter"),
77  multiple=True,
78  )
79  dcrNImages = pipeBase.connectionTypes.Output(
80  doc="Output image of number of input images per pixel",
81  name="{outputCoaddName}Coadd_nImage",
82  storageClass="ImageU",
83  dimensions=("tract", "patch", "skymap", "band", "subfilter"),
84  multiple=True,
85  )
86 
87  def __init__(self, *, config=None):
88  super().__init__(config=config)
89  if not config.doWrite:
90  self.outputs.remove("dcrCoadds")
91  if not config.doNImage:
92  self.outputs.remove("dcrNImages")
93  # Remove outputs inherited from ``AssembleCoaddConnections`` that are not used
94  self.outputs.remove("coaddExposure")
95  self.outputs.remove("nImage")
96 
97 
98 class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig,
99  pipelineConnections=DcrAssembleCoaddConnections):
100  dcrNumSubfilters = pexConfig.Field(
101  dtype=int,
102  doc="Number of sub-filters to forward model chromatic effects to fit the supplied exposures.",
103  default=3,
104  )
105  maxNumIter = pexConfig.Field(
106  dtype=int,
107  optional=True,
108  doc="Maximum number of iterations of forward modeling.",
109  default=None,
110  )
111  minNumIter = pexConfig.Field(
112  dtype=int,
113  optional=True,
114  doc="Minimum number of iterations of forward modeling.",
115  default=None,
116  )
117  convergenceThreshold = pexConfig.Field(
118  dtype=float,
119  doc="Target relative change in convergence between iterations of forward modeling.",
120  default=0.001,
121  )
122  useConvergence = pexConfig.Field(
123  dtype=bool,
124  doc="Use convergence test as a forward modeling end condition?"
125  "If not set, skips calculating convergence and runs for ``maxNumIter`` iterations",
126  default=True,
127  )
128  baseGain = pexConfig.Field(
129  dtype=float,
130  optional=True,
131  doc="Relative weight to give the new solution vs. the last solution when updating the model."
132  "A value of 1.0 gives equal weight to both solutions."
133  "Small values imply slower convergence of the solution, but can "
134  "help prevent overshooting and failures in the fit."
135  "If ``baseGain`` is None, a conservative gain "
136  "will be calculated from the number of subfilters. ",
137  default=None,
138  )
139  useProgressiveGain = pexConfig.Field(
140  dtype=bool,
141  doc="Use a gain that slowly increases above ``baseGain`` to accelerate convergence? "
142  "When calculating the next gain, we use up to 5 previous gains and convergence values."
143  "Can be set to False to force the model to change at the rate of ``baseGain``. ",
144  default=True,
145  )
146  doAirmassWeight = pexConfig.Field(
147  dtype=bool,
148  doc="Weight exposures by airmass? Useful if there are relatively few high-airmass observations.",
149  default=False,
150  )
151  modelWeightsWidth = pexConfig.Field(
152  dtype=float,
153  doc="Width of the region around detected sources to include in the DcrModel.",
154  default=3,
155  )
156  useModelWeights = pexConfig.Field(
157  dtype=bool,
158  doc="Width of the region around detected sources to include in the DcrModel.",
159  default=True,
160  )
161  splitSubfilters = pexConfig.Field(
162  dtype=bool,
163  doc="Calculate DCR for two evenly-spaced wavelengths in each subfilter."
164  "Instead of at the midpoint",
165  default=True,
166  )
167  splitThreshold = pexConfig.Field(
168  dtype=float,
169  doc="Minimum DCR difference within a subfilter to use ``splitSubfilters``, in pixels."
170  "Set to 0 to always split the subfilters.",
171  default=0.1,
172  )
173  regularizeModelIterations = pexConfig.Field(
174  dtype=float,
175  doc="Maximum relative change of the model allowed between iterations."
176  "Set to zero to disable.",
177  default=2.,
178  )
179  regularizeModelFrequency = pexConfig.Field(
180  dtype=float,
181  doc="Maximum relative change of the model allowed between subfilters."
182  "Set to zero to disable.",
183  default=4.,
184  )
185  convergenceMaskPlanes = pexConfig.ListField(
186  dtype=str,
187  default=["DETECTED"],
188  doc="Mask planes to use to calculate convergence."
189  )
190  regularizationWidth = pexConfig.Field(
191  dtype=int,
192  default=2,
193  doc="Minimum radius of a region to include in regularization, in pixels."
194  )
195  imageInterpOrder = pexConfig.Field(
196  dtype=int,
197  doc="The order of the spline interpolation used to shift the image plane.",
198  default=3,
199  )
200  accelerateModel = pexConfig.Field(
201  dtype=float,
202  doc="Factor to amplify the differences between model planes by to speed convergence.",
203  default=3,
204  )
205  doCalculatePsf = pexConfig.Field(
206  dtype=bool,
207  doc="Set to detect stars and recalculate the PSF from the final coadd."
208  "Otherwise the PSF is estimated from a selection of the best input exposures",
209  default=False,
210  )
211  detectPsfSources = pexConfig.ConfigurableField(
212  target=measAlg.SourceDetectionTask,
213  doc="Task to detect sources for PSF measurement, if ``doCalculatePsf`` is set.",
214  )
215  measurePsfSources = pexConfig.ConfigurableField(
216  target=SingleFrameMeasurementTask,
217  doc="Task to measure sources for PSF measurement, if ``doCalculatePsf`` is set."
218  )
219  measurePsf = pexConfig.ConfigurableField(
220  target=MeasurePsfTask,
221  doc="Task to measure the PSF of the coadd, if ``doCalculatePsf`` is set.",
222  )
223  effectiveWavelength = pexConfig.Field(
224  doc="Effective wavelength of the filter, in nm."
225  "Required if transmission curves aren't used."
226  "Support for using transmission curves is to be added in DM-13668.",
227  dtype=float,
228  )
229  bandwidth = pexConfig.Field(
230  doc="Bandwidth of the physical filter, in nm."
231  "Required if transmission curves aren't used."
232  "Support for using transmission curves is to be added in DM-13668.",
233  dtype=float,
234  )
235 
236  def setDefaults(self):
237  CompareWarpAssembleCoaddConfig.setDefaults(self)
238  self.assembleStaticSkyModel.retarget(CompareWarpAssembleCoaddTask)
239  self.doNImage = True
240  self.assembleStaticSkyModel.warpType = self.warpType
241  # The deepCoadd and nImage files will be overwritten by this Task, so don't write them the first time
242  self.assembleStaticSkyModel.doNImage = False
243  self.assembleStaticSkyModel.doWrite = False
244  self.detectPsfSources.returnOriginalFootprints = False
245  self.detectPsfSources.thresholdPolarity = "positive"
246  # Only use bright sources for PSF measurement
247  self.detectPsfSources.thresholdValue = 50
248  self.detectPsfSources.nSigmaToGrow = 2
249  # A valid star for PSF measurement should at least fill 5x5 pixels
250  self.detectPsfSources.minPixels = 25
251  # Use the variance plane to calculate signal to noise
252  self.detectPsfSources.thresholdType = "pixel_stdev"
253  # The signal to noise limit is good enough, while the flux limit is set
254  # in dimensionless units and may not be appropriate for all data sets.
255  self.measurePsf.starSelector["objectSize"].doFluxLimit = False
256 
257 
258 class DcrAssembleCoaddTask(CompareWarpAssembleCoaddTask):
259  """Assemble DCR coadded images from a set of warps.
260 
261  Attributes
262  ----------
263  bufferSize : `int`
264  The number of pixels to grow each subregion by to allow for DCR.
265 
266  Notes
267  -----
268  As with AssembleCoaddTask, we want to assemble a coadded image from a set of
269  Warps (also called coadded temporary exposures), including the effects of
270  Differential Chromatic Refraction (DCR).
271  For full details of the mathematics and algorithm, please see
272  DMTN-037: DCR-matched template generation (https://dmtn-037.lsst.io).
273 
274  This Task produces a DCR-corrected deepCoadd, as well as a dcrCoadd for
275  each subfilter used in the iterative calculation.
276  It begins by dividing the bandpass-defining filter into N equal bandwidth
277  "subfilters", and divides the flux in each pixel from an initial coadd
278  equally into each as a "dcrModel". Because the airmass and parallactic
279  angle of each individual exposure is known, we can calculate the shift
280  relative to the center of the band in each subfilter due to DCR. For each
281  exposure we apply this shift as a linear transformation to the dcrModels
282  and stack the results to produce a DCR-matched exposure. The matched
283  exposures are subtracted from the input exposures to produce a set of
284  residual images, and these residuals are reverse shifted for each
285  exposures' subfilters and stacked. The shifted and stacked residuals are
286  added to the dcrModels to produce a new estimate of the flux in each pixel
287  within each subfilter. The dcrModels are solved for iteratively, which
288  continues until the solution from a new iteration improves by less than
289  a set percentage, or a maximum number of iterations is reached.
290  Two forms of regularization are employed to reduce unphysical results.
291  First, the new solution is averaged with the solution from the previous
292  iteration, which mitigates oscillating solutions where the model
293  overshoots with alternating very high and low values.
294  Second, a common degeneracy when the data have a limited range of airmass or
295  parallactic angle values is for one subfilter to be fit with very low or
296  negative values, while another subfilter is fit with very high values. This
297  typically appears in the form of holes next to sources in one subfilter,
298  and corresponding extended wings in another. Because each subfilter has
299  a narrow bandwidth we assume that physical sources that are above the noise
300  level will not vary in flux by more than a factor of `frequencyClampFactor`
301  between subfilters, and pixels that have flux deviations larger than that
302  factor will have the excess flux distributed evenly among all subfilters.
303  If `splitSubfilters` is set, then each subfilter will be further sub-
304  divided during the forward modeling step (only). This approximates using
305  a higher number of subfilters that may be necessary for high airmass
306  observations, but does not increase the number of free parameters in the
307  fit. This is needed when there are high airmass observations which would
308  otherwise have significant DCR even within a subfilter. Because calculating
309  the shifted images takes most of the time, splitting the subfilters is
310  turned off by way of the `splitThreshold` option for low-airmass
311  observations that do not suffer from DCR within a subfilter.
312  """
313 
314  ConfigClass = DcrAssembleCoaddConfig
315  _DefaultName = "dcrAssembleCoadd"
316 
317  def __init__(self, *args, **kwargs):
318  super().__init__(*args, **kwargs)
319  if self.config.doCalculatePsf:
320  self.schema = afwTable.SourceTable.makeMinimalSchema()
321  self.makeSubtask("detectPsfSources", schema=self.schema)
322  self.makeSubtask("measurePsfSources", schema=self.schema)
323  self.makeSubtask("measurePsf", schema=self.schema)
324 
325  @utils.inheritDoc(pipeBase.PipelineTask)
326  def runQuantum(self, butlerQC, inputRefs, outputRefs):
327  # Docstring to be formatted with info from PipelineTask.runQuantum
328  """
329  Notes
330  -----
331  Assemble a coadd from a set of Warps.
332 
333  PipelineTask (Gen3) entry point to Coadd a set of Warps.
334  Analogous to `runDataRef`, it prepares all the data products to be
335  passed to `run`, and processes the results before returning a struct
336  of results to be written out. AssembleCoadd cannot fit all Warps in memory.
337  Therefore, its inputs are accessed subregion by subregion
338  by the Gen3 `DeferredDatasetHandle` that is analagous to the Gen2
339  `lsst.daf.persistence.ButlerDataRef`. Any updates to this method should
340  correspond to an update in `runDataRef` while both entry points
341  are used.
342  """
343  inputData = butlerQC.get(inputRefs)
344 
345  # Construct skyInfo expected by run
346  # Do not remove skyMap from inputData in case makeSupplementaryDataGen3 needs it
347  skyMap = inputData["skyMap"]
348  outputDataId = butlerQC.quantum.dataId
349 
350  inputData['skyInfo'] = makeSkyInfo(skyMap,
351  tractId=outputDataId['tract'],
352  patchId=outputDataId['patch'])
353 
354  # Construct list of input Deferred Datasets
355  # These quack a bit like like Gen2 DataRefs
356  warpRefList = inputData['inputWarps']
357  # Perform same middle steps as `runDataRef` does
358  inputs = self.prepareInputs(warpRefList)
359  self.log.info("Found %d %s", len(inputs.tempExpRefList),
360  self.getTempExpDatasetName(self.warpType))
361  if len(inputs.tempExpRefList) == 0:
362  self.log.warning("No coadd temporary exposures found")
363  return
364 
365  supplementaryData = self.makeSupplementaryDataGen3(butlerQC, inputRefs, outputRefs)
366  retStruct = self.run(inputData['skyInfo'], inputs.tempExpRefList, inputs.imageScalerList,
367  inputs.weightList, supplementaryData=supplementaryData)
368 
369  inputData.setdefault('brightObjectMask', None)
370  for subfilter in range(self.config.dcrNumSubfilters):
371  # Use the PSF of the stacked dcrModel, and do not recalculate the PSF for each subfilter
372  retStruct.dcrCoadds[subfilter].setPsf(retStruct.coaddExposure.getPsf())
373  self.processResults(retStruct.dcrCoadds[subfilter], inputData['brightObjectMask'], outputDataId)
374 
375  if self.config.doWrite:
376  butlerQC.put(retStruct, outputRefs)
377  return retStruct
378 
379  @timeMethod
380  def runDataRef(self, dataRef, selectDataList=None, warpRefList=None):
381  """Assemble a coadd from a set of warps.
382 
383  Coadd a set of Warps. Compute weights to be applied to each Warp and
384  find scalings to match the photometric zeropoint to a reference Warp.
385  Assemble the Warps using run method.
386  Forward model chromatic effects across multiple subfilters,
387  and subtract from the input Warps to build sets of residuals.
388  Use the residuals to construct a new ``DcrModel`` for each subfilter,
389  and iterate until the model converges.
390  Interpolate over NaNs and optionally write the coadd to disk.
391  Return the coadded exposure.
392 
393  Parameters
394  ----------
395  dataRef : `lsst.daf.persistence.ButlerDataRef`
396  Data reference defining the patch for coaddition and the
397  reference Warp
398  selectDataList : `list` of `lsst.daf.persistence.ButlerDataRef`
399  List of data references to warps. Data to be coadded will be
400  selected from this list based on overlap with the patch defined by
401  the data reference.
402 
403  Returns
404  -------
405  results : `lsst.pipe.base.Struct`
406  The Struct contains the following fields:
407 
408  - ``coaddExposure``: coadded exposure (`lsst.afw.image.Exposure`)
409  - ``nImage``: exposure count image (`lsst.afw.image.ImageU`)
410  - ``dcrCoadds``: `list` of coadded exposures for each subfilter
411  - ``dcrNImages``: `list` of exposure count images for each subfilter
412  """
413  if (selectDataList is None and warpRefList is None) or (selectDataList and warpRefList):
414  raise RuntimeError("runDataRef must be supplied either a selectDataList or warpRefList")
415 
416  skyInfo = self.getSkyInfo(dataRef)
417  if warpRefList is None:
418  calExpRefList = self.selectExposures(dataRef, skyInfo, selectDataList=selectDataList)
419  if len(calExpRefList) == 0:
420  self.log.warning("No exposures to coadd")
421  return
422  self.log.info("Coadding %d exposures", len(calExpRefList))
423 
424  warpRefList = self.getTempExpRefList(dataRef, calExpRefList)
425 
426  inputData = self.prepareInputs(warpRefList)
427  self.log.info("Found %d %s", len(inputData.tempExpRefList),
428  self.getTempExpDatasetName(self.warpType))
429  if len(inputData.tempExpRefList) == 0:
430  self.log.warning("No coadd temporary exposures found")
431  return
432 
433  supplementaryData = self.makeSupplementaryData(dataRef, warpRefList=inputData.tempExpRefList)
434 
435  results = self.run(skyInfo, inputData.tempExpRefList, inputData.imageScalerList,
436  inputData.weightList, supplementaryData=supplementaryData)
437  if results is None:
438  self.log.warning("Could not construct DcrModel for patch %s: no data to coadd.",
439  skyInfo.patchInfo.getIndex())
440  return
441 
442  if self.config.doCalculatePsf:
443  self.measureCoaddPsf(results.coaddExposure)
444  brightObjects = self.readBrightObjectMasks(dataRef) if self.config.doMaskBrightObjects else None
445  for subfilter in range(self.config.dcrNumSubfilters):
446  # Use the PSF of the stacked dcrModel, and do not recalculate the PSF for each subfilter
447  results.dcrCoadds[subfilter].setPsf(results.coaddExposure.getPsf())
448  self.processResults(results.dcrCoadds[subfilter],
449  brightObjectMasks=brightObjects, dataId=dataRef.dataId)
450  if self.config.doWrite:
451  self.log.info("Persisting dcrCoadd")
452  dataRef.put(results.dcrCoadds[subfilter], "dcrCoadd", subfilter=subfilter,
453  numSubfilters=self.config.dcrNumSubfilters)
454  if self.config.doNImage and results.dcrNImages is not None:
455  dataRef.put(results.dcrNImages[subfilter], "dcrCoadd_nImage", subfilter=subfilter,
456  numSubfilters=self.config.dcrNumSubfilters)
457 
458  return results
459 
460  @utils.inheritDoc(AssembleCoaddTask)
461  def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs):
462  """Load the previously-generated template coadd.
463 
464  This can be removed entirely once we no longer support the Gen 2 butler.
465 
466  Returns
467  -------
468  templateCoadd : `lsst.pipe.base.Struct`
469  Result struct with components:
470 
471  - ``templateCoadd``: coadded exposure (`lsst.afw.image.ExposureF`)
472  """
473  templateCoadd = butlerQC.get(inputRefs.templateExposure)
474 
475  return pipeBase.Struct(templateCoadd=templateCoadd)
476 
477  def measureCoaddPsf(self, coaddExposure):
478  """Detect sources on the coadd exposure and measure the final PSF.
479 
480  Parameters
481  ----------
482  coaddExposure : `lsst.afw.image.Exposure`
483  The final coadded exposure.
484  """
485  table = afwTable.SourceTable.make(self.schema)
486  detResults = self.detectPsfSources.run(table, coaddExposure, clearMask=False)
487  coaddSources = detResults.sources
488  self.measurePsfSources.run(
489  measCat=coaddSources,
490  exposure=coaddExposure
491  )
492  # Measure the PSF on the stacked subfilter coadds if possible.
493  # We should already have a decent estimate of the coadd PSF, however,
494  # so in case of any errors simply log them as a warning and use the
495  # default PSF.
496  try:
497  psfResults = self.measurePsf.run(coaddExposure, coaddSources)
498  except Exception as e:
499  self.log.warning("Unable to calculate PSF, using default coadd PSF: %s", e)
500  else:
501  coaddExposure.setPsf(psfResults.psf)
502 
503  def prepareDcrInputs(self, templateCoadd, warpRefList, weightList):
504  """Prepare the DCR coadd by iterating through the visitInfo of the input warps.
505 
506  Sets the property ``bufferSize``.
507 
508  Parameters
509  ----------
510  templateCoadd : `lsst.afw.image.ExposureF`
511  The initial coadd exposure before accounting for DCR.
512  warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
513  `lsst.daf.persistence.ButlerDataRef`
514  The data references to the input warped exposures.
515  weightList : `list` of `float`
516  The weight to give each input exposure in the coadd
517  Will be modified in place if ``doAirmassWeight`` is set.
518 
519  Returns
520  -------
521  dcrModels : `lsst.pipe.tasks.DcrModel`
522  Best fit model of the true sky after correcting chromatic effects.
523 
524  Raises
525  ------
526  NotImplementedError
527  If ``lambdaMin`` is missing from the Mapper class of the obs package being used.
528  """
529  sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
530  filterLabel = templateCoadd.getFilterLabel()
531  tempExpName = self.getTempExpDatasetName(self.warpType)
532  dcrShifts = []
533  airmassDict = {}
534  angleDict = {}
535  psfSizeDict = {}
536  for visitNum, warpExpRef in enumerate(warpRefList):
537  if isinstance(warpExpRef, DeferredDatasetHandle):
538  # Gen 3 API
539  visitInfo = warpExpRef.get(component="visitInfo")
540  psf = warpExpRef.get(component="psf")
541  else:
542  # Gen 2 API. Delete this when Gen 2 retired
543  visitInfo = warpExpRef.get(tempExpName + "_visitInfo")
544  psf = warpExpRef.get(tempExpName).getPsf()
545  visit = warpExpRef.dataId["visit"]
546  psfSize = psf.computeShape().getDeterminantRadius()*sigma2fwhm
547  airmass = visitInfo.getBoresightAirmass()
548  parallacticAngle = visitInfo.getBoresightParAngle().asDegrees()
549  airmassDict[visit] = airmass
550  angleDict[visit] = parallacticAngle
551  psfSizeDict[visit] = psfSize
552  if self.config.doAirmassWeight:
553  weightList[visitNum] *= airmass
554  dcrShifts.append(np.max(np.abs(calculateDcr(visitInfo, templateCoadd.getWcs(),
555  self.config.effectiveWavelength,
556  self.config.bandwidth,
557  self.config.dcrNumSubfilters))))
558  self.log.info("Selected airmasses:\n%s", airmassDict)
559  self.log.info("Selected parallactic angles:\n%s", angleDict)
560  self.log.info("Selected PSF sizes:\n%s", psfSizeDict)
561  self.bufferSize = int(np.ceil(np.max(dcrShifts)) + 1)
562  try:
563  psf = self.selectCoaddPsf(templateCoadd, warpRefList)
564  except Exception as e:
565  self.log.warning("Unable to calculate restricted PSF, using default coadd PSF: %s", e)
566  else:
567  psf = templateCoadd.getPsf()
568  dcrModels = DcrModel.fromImage(templateCoadd.maskedImage,
569  self.config.dcrNumSubfilters,
570  effectiveWavelength=self.config.effectiveWavelength,
571  bandwidth=self.config.bandwidth,
572  filterLabel=filterLabel,
573  psf=psf)
574  return dcrModels
575 
576  @timeMethod
577  def run(self, skyInfo, warpRefList, imageScalerList, weightList,
578  supplementaryData=None):
579  """Assemble the coadd.
580 
581  Requires additional inputs Struct ``supplementaryData`` to contain a
582  ``templateCoadd`` that serves as the model of the static sky.
583 
584  Find artifacts and apply them to the warps' masks creating a list of
585  alternative masks with a new "CLIPPED" plane and updated "NO_DATA" plane
586  Then pass these alternative masks to the base class's assemble method.
587 
588  Divide the ``templateCoadd`` evenly between each subfilter of a
589  ``DcrModel`` as the starting best estimate of the true wavelength-
590  dependent sky. Forward model the ``DcrModel`` using the known
591  chromatic effects in each subfilter and calculate a convergence metric
592  based on how well the modeled template matches the input warps. If
593  the convergence has not yet reached the desired threshold, then shift
594  and stack the residual images to build a new ``DcrModel``. Apply
595  conditioning to prevent oscillating solutions between iterations or
596  between subfilters.
597 
598  Once the ``DcrModel`` reaches convergence or the maximum number of
599  iterations has been reached, fill the metadata for each subfilter
600  image and make them proper ``coaddExposure``s.
601 
602  Parameters
603  ----------
604  skyInfo : `lsst.pipe.base.Struct`
605  Patch geometry information, from getSkyInfo
606  warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
607  `lsst.daf.persistence.ButlerDataRef`
608  The data references to the input warped exposures.
609  imageScalerList : `list` of `lsst.pipe.task.ImageScaler`
610  The image scalars correct for the zero point of the exposures.
611  weightList : `list` of `float`
612  The weight to give each input exposure in the coadd
613  supplementaryData : `lsst.pipe.base.Struct`
614  Result struct returned by ``makeSupplementaryData`` with components:
615 
616  - ``templateCoadd``: coadded exposure (`lsst.afw.image.Exposure`)
617 
618  Returns
619  -------
620  result : `lsst.pipe.base.Struct`
621  Result struct with components:
622 
623  - ``coaddExposure``: coadded exposure (`lsst.afw.image.Exposure`)
624  - ``nImage``: exposure count image (`lsst.afw.image.ImageU`)
625  - ``dcrCoadds``: `list` of coadded exposures for each subfilter
626  - ``dcrNImages``: `list` of exposure count images for each subfilter
627  """
628  minNumIter = self.config.minNumIter or self.config.dcrNumSubfilters
629  maxNumIter = self.config.maxNumIter or self.config.dcrNumSubfilters*3
630  templateCoadd = supplementaryData.templateCoadd
631  baseMask = templateCoadd.mask.clone()
632  # The variance plane is for each subfilter
633  # and should be proportionately lower than the full-band image
634  baseVariance = templateCoadd.variance.clone()
635  baseVariance /= self.config.dcrNumSubfilters
636  spanSetMaskList = self.findArtifacts(templateCoadd, warpRefList, imageScalerList)
637  # Note that the mask gets cleared in ``findArtifacts``, but we want to preserve the mask.
638  templateCoadd.setMask(baseMask)
639  badMaskPlanes = self.config.badMaskPlanes[:]
640  # Note that is important that we do not add "CLIPPED" to ``badMaskPlanes``
641  # This is because pixels in observations that are significantly affect by DCR
642  # are likely to have many pixels that are both "DETECTED" and "CLIPPED",
643  # but those are necessary to constrain the DCR model.
644  badPixelMask = templateCoadd.mask.getPlaneBitMask(badMaskPlanes)
645 
646  stats = self.prepareStats(mask=badPixelMask)
647  dcrModels = self.prepareDcrInputs(templateCoadd, warpRefList, weightList)
648  if self.config.doNImage:
649  dcrNImages, dcrWeights = self.calculateNImage(dcrModels, skyInfo.bbox, warpRefList,
650  spanSetMaskList, stats.ctrl)
651  nImage = afwImage.ImageU(skyInfo.bbox)
652  # Note that this nImage will be a factor of dcrNumSubfilters higher than
653  # the nImage returned by assembleCoadd for most pixels. This is because each
654  # subfilter may have a different nImage, and fractional values are not allowed.
655  for dcrNImage in dcrNImages:
656  nImage += dcrNImage
657  else:
658  dcrNImages = None
659 
660  subregionSize = geom.Extent2I(*self.config.subregionSize)
661  nSubregions = (ceil(skyInfo.bbox.getHeight()/subregionSize[1])
662  * ceil(skyInfo.bbox.getWidth()/subregionSize[0]))
663  subIter = 0
664  for subBBox in self._subBBoxIter(skyInfo.bbox, subregionSize):
665  modelIter = 0
666  subIter += 1
667  self.log.info("Computing coadd over patch %s subregion %s of %s: %s",
668  skyInfo.patchInfo.getIndex(), subIter, nSubregions, subBBox)
669  dcrBBox = geom.Box2I(subBBox)
670  dcrBBox.grow(self.bufferSize)
671  dcrBBox.clip(dcrModels.bbox)
672  modelWeights = self.calculateModelWeights(dcrModels, dcrBBox)
673  subExposures = self.loadSubExposures(dcrBBox, stats.ctrl, warpRefList,
674  imageScalerList, spanSetMaskList)
675  convergenceMetric = self.calculateConvergence(dcrModels, subExposures, subBBox,
676  warpRefList, weightList, stats.ctrl)
677  self.log.info("Initial convergence : %s", convergenceMetric)
678  convergenceList = [convergenceMetric]
679  gainList = []
680  convergenceCheck = 1.
681  refImage = templateCoadd.image
682  while (convergenceCheck > self.config.convergenceThreshold or modelIter <= minNumIter):
683  gain = self.calculateGain(convergenceList, gainList)
684  self.dcrAssembleSubregion(dcrModels, subExposures, subBBox, dcrBBox, warpRefList,
685  stats.ctrl, convergenceMetric, gain,
686  modelWeights, refImage, dcrWeights)
687  if self.config.useConvergence:
688  convergenceMetric = self.calculateConvergence(dcrModels, subExposures, subBBox,
689  warpRefList, weightList, stats.ctrl)
690  if convergenceMetric == 0:
691  self.log.warning("Coadd patch %s subregion %s had convergence metric of 0.0 which is "
692  "most likely due to there being no valid data in the region.",
693  skyInfo.patchInfo.getIndex(), subIter)
694  break
695  convergenceCheck = (convergenceList[-1] - convergenceMetric)/convergenceMetric
696  if (convergenceCheck < 0) & (modelIter > minNumIter):
697  self.log.warning("Coadd patch %s subregion %s diverged before reaching maximum "
698  "iterations or desired convergence improvement of %s."
699  " Divergence: %s",
700  skyInfo.patchInfo.getIndex(), subIter,
701  self.config.convergenceThreshold, convergenceCheck)
702  break
703  convergenceList.append(convergenceMetric)
704  if modelIter > maxNumIter:
705  if self.config.useConvergence:
706  self.log.warning("Coadd patch %s subregion %s reached maximum iterations "
707  "before reaching desired convergence improvement of %s."
708  " Final convergence improvement: %s",
709  skyInfo.patchInfo.getIndex(), subIter,
710  self.config.convergenceThreshold, convergenceCheck)
711  break
712 
713  if self.config.useConvergence:
714  self.log.info("Iteration %s with convergence metric %s, %.4f%% improvement (gain: %.2f)",
715  modelIter, convergenceMetric, 100.*convergenceCheck, gain)
716  modelIter += 1
717  else:
718  if self.config.useConvergence:
719  self.log.info("Coadd patch %s subregion %s finished with "
720  "convergence metric %s after %s iterations",
721  skyInfo.patchInfo.getIndex(), subIter, convergenceMetric, modelIter)
722  else:
723  self.log.info("Coadd patch %s subregion %s finished after %s iterations",
724  skyInfo.patchInfo.getIndex(), subIter, modelIter)
725  if self.config.useConvergence and convergenceMetric > 0:
726  self.log.info("Final convergence improvement was %.4f%% overall",
727  100*(convergenceList[0] - convergenceMetric)/convergenceMetric)
728 
729  dcrCoadds = self.fillCoadd(dcrModels, skyInfo, warpRefList, weightList,
730  calibration=self.scaleZeroPoint.getPhotoCalib(),
731  coaddInputs=templateCoadd.getInfo().getCoaddInputs(),
732  mask=baseMask,
733  variance=baseVariance)
734  coaddExposure = self.stackCoadd(dcrCoadds)
735  return pipeBase.Struct(coaddExposure=coaddExposure, nImage=nImage,
736  dcrCoadds=dcrCoadds, dcrNImages=dcrNImages)
737 
738  def calculateNImage(self, dcrModels, bbox, warpRefList, spanSetMaskList, statsCtrl):
739  """Calculate the number of exposures contributing to each subfilter.
740 
741  Parameters
742  ----------
743  dcrModels : `lsst.pipe.tasks.DcrModel`
744  Best fit model of the true sky after correcting chromatic effects.
745  bbox : `lsst.geom.box.Box2I`
746  Bounding box of the patch to coadd.
747  warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
748  `lsst.daf.persistence.ButlerDataRef`
749  The data references to the input warped exposures.
750  spanSetMaskList : `list` of `dict` containing spanSet lists, or None
751  Each element of the `dict` contains the new mask plane name
752  (e.g. "CLIPPED and/or "NO_DATA") as the key,
753  and the list of SpanSets to apply to the mask.
754  statsCtrl : `lsst.afw.math.StatisticsControl`
755  Statistics control object for coadd
756 
757  Returns
758  -------
759  dcrNImages : `list` of `lsst.afw.image.ImageU`
760  List of exposure count images for each subfilter
761  dcrWeights : `list` of `lsst.afw.image.ImageF`
762  Per-pixel weights for each subfilter.
763  Equal to 1/(number of unmasked images contributing to each pixel).
764  """
765  dcrNImages = [afwImage.ImageU(bbox) for subfilter in range(self.config.dcrNumSubfilters)]
766  dcrWeights = [afwImage.ImageF(bbox) for subfilter in range(self.config.dcrNumSubfilters)]
767  tempExpName = self.getTempExpDatasetName(self.warpType)
768  for warpExpRef, altMaskSpans in zip(warpRefList, spanSetMaskList):
769  if isinstance(warpExpRef, DeferredDatasetHandle):
770  # Gen 3 API
771  exposure = warpExpRef.get(parameters={'bbox': bbox})
772  else:
773  # Gen 2 API. Delete this when Gen 2 retired
774  exposure = warpExpRef.get(tempExpName + "_sub", bbox=bbox)
775  visitInfo = exposure.getInfo().getVisitInfo()
776  wcs = exposure.getInfo().getWcs()
777  mask = exposure.mask
778  if altMaskSpans is not None:
779  self.applyAltMaskPlanes(mask, altMaskSpans)
780  weightImage = np.zeros_like(exposure.image.array)
781  weightImage[(mask.array & statsCtrl.getAndMask()) == 0] = 1.
782  # The weights must be shifted in exactly the same way as the residuals,
783  # because they will be used as the denominator in the weighted average of residuals.
784  weightsGenerator = self.dcrResiduals(weightImage, visitInfo, wcs,
785  dcrModels.effectiveWavelength, dcrModels.bandwidth)
786  for shiftedWeights, dcrNImage, dcrWeight in zip(weightsGenerator, dcrNImages, dcrWeights):
787  dcrNImage.array += np.rint(shiftedWeights).astype(dcrNImage.array.dtype)
788  dcrWeight.array += shiftedWeights
789  # Exclude any pixels that don't have at least one exposure contributing in all subfilters
790  weightsThreshold = 1.
791  goodPix = dcrWeights[0].array > weightsThreshold
792  for weights in dcrWeights[1:]:
793  goodPix = (weights.array > weightsThreshold) & goodPix
794  for subfilter in range(self.config.dcrNumSubfilters):
795  dcrWeights[subfilter].array[goodPix] = 1./dcrWeights[subfilter].array[goodPix]
796  dcrWeights[subfilter].array[~goodPix] = 0.
797  dcrNImages[subfilter].array[~goodPix] = 0
798  return (dcrNImages, dcrWeights)
799 
800  def dcrAssembleSubregion(self, dcrModels, subExposures, bbox, dcrBBox, warpRefList,
801  statsCtrl, convergenceMetric,
802  gain, modelWeights, refImage, dcrWeights):
803  """Assemble the DCR coadd for a sub-region.
804 
805  Build a DCR-matched template for each input exposure, then shift the
806  residuals according to the DCR in each subfilter.
807  Stack the shifted residuals and apply them as a correction to the
808  solution from the previous iteration.
809  Restrict the new model solutions from varying by more than a factor of
810  `modelClampFactor` from the last solution, and additionally restrict the
811  individual subfilter models from varying by more than a factor of
812  `frequencyClampFactor` from their average.
813  Finally, mitigate potentially oscillating solutions by averaging the new
814  solution with the solution from the previous iteration, weighted by
815  their convergence metric.
816 
817  Parameters
818  ----------
819  dcrModels : `lsst.pipe.tasks.DcrModel`
820  Best fit model of the true sky after correcting chromatic effects.
821  subExposures : `dict` of `lsst.afw.image.ExposureF`
822  The pre-loaded exposures for the current subregion.
823  bbox : `lsst.geom.box.Box2I`
824  Bounding box of the subregion to coadd.
825  dcrBBox : `lsst.geom.box.Box2I`
826  Sub-region of the coadd which includes a buffer to allow for DCR.
827  warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
828  `lsst.daf.persistence.ButlerDataRef`
829  The data references to the input warped exposures.
830  statsCtrl : `lsst.afw.math.StatisticsControl`
831  Statistics control object for coadd
832  convergenceMetric : `float`
833  Quality of fit metric for the matched templates of the input images.
834  gain : `float`, optional
835  Relative weight to give the new solution when updating the model.
836  modelWeights : `numpy.ndarray` or `float`
837  A 2D array of weight values that tapers smoothly to zero away from detected sources.
838  Set to a placeholder value of 1.0 if ``self.config.useModelWeights`` is False.
839  refImage : `lsst.afw.image.Image`
840  A reference image used to supply the default pixel values.
841  dcrWeights : `list` of `lsst.afw.image.Image`
842  Per-pixel weights for each subfilter.
843  Equal to 1/(number of unmasked images contributing to each pixel).
844  """
845  residualGeneratorList = []
846 
847  for warpExpRef in warpRefList:
848  visit = warpExpRef.dataId["visit"]
849  exposure = subExposures[visit]
850  visitInfo = exposure.getInfo().getVisitInfo()
851  wcs = exposure.getInfo().getWcs()
852  templateImage = dcrModels.buildMatchedTemplate(exposure=exposure,
853  order=self.config.imageInterpOrder,
854  splitSubfilters=self.config.splitSubfilters,
855  splitThreshold=self.config.splitThreshold,
856  amplifyModel=self.config.accelerateModel)
857  residual = exposure.image.array - templateImage.array
858  # Note that the variance plane here is used to store weights, not the actual variance
859  residual *= exposure.variance.array
860  # The residuals are stored as a list of generators.
861  # This allows the residual for a given subfilter and exposure to be created
862  # on the fly, instead of needing to store them all in memory.
863  residualGeneratorList.append(self.dcrResiduals(residual, visitInfo, wcs,
864  dcrModels.effectiveWavelength,
865  dcrModels.bandwidth))
866 
867  dcrSubModelOut = self.newModelFromResidual(dcrModels, residualGeneratorList, dcrBBox, statsCtrl,
868  gain=gain,
869  modelWeights=modelWeights,
870  refImage=refImage,
871  dcrWeights=dcrWeights)
872  dcrModels.assign(dcrSubModelOut, bbox)
873 
874  def dcrResiduals(self, residual, visitInfo, wcs, effectiveWavelength, bandwidth):
875  """Prepare a residual image for stacking in each subfilter by applying the reverse DCR shifts.
876 
877  Parameters
878  ----------
879  residual : `numpy.ndarray`
880  The residual masked image for one exposure,
881  after subtracting the matched template
882  visitInfo : `lsst.afw.image.VisitInfo`
883  Metadata for the exposure.
884  wcs : `lsst.afw.geom.SkyWcs`
885  Coordinate system definition (wcs) for the exposure.
886 
887  Yields
888  ------
889  residualImage : `numpy.ndarray`
890  The residual image for the next subfilter, shifted for DCR.
891  """
892  # Pre-calculate the spline-filtered residual image, so that step can be
893  # skipped in the shift calculation in `applyDcr`.
894  filteredResidual = ndimage.spline_filter(residual, order=self.config.imageInterpOrder)
895  # Note that `splitSubfilters` is always turned off in the reverse direction.
896  # This option introduces additional blurring if applied to the residuals.
897  dcrShift = calculateDcr(visitInfo, wcs, effectiveWavelength, bandwidth, self.config.dcrNumSubfilters,
898  splitSubfilters=False)
899  for dcr in dcrShift:
900  yield applyDcr(filteredResidual, dcr, useInverse=True, splitSubfilters=False,
901  doPrefilter=False, order=self.config.imageInterpOrder)
902 
903  def newModelFromResidual(self, dcrModels, residualGeneratorList, dcrBBox, statsCtrl,
904  gain, modelWeights, refImage, dcrWeights):
905  """Calculate a new DcrModel from a set of image residuals.
906 
907  Parameters
908  ----------
909  dcrModels : `lsst.pipe.tasks.DcrModel`
910  Current model of the true sky after correcting chromatic effects.
911  residualGeneratorList : `generator` of `numpy.ndarray`
912  The residual image for the next subfilter, shifted for DCR.
913  dcrBBox : `lsst.geom.box.Box2I`
914  Sub-region of the coadd which includes a buffer to allow for DCR.
915  statsCtrl : `lsst.afw.math.StatisticsControl`
916  Statistics control object for coadd
917  gain : `float`
918  Relative weight to give the new solution when updating the model.
919  modelWeights : `numpy.ndarray` or `float`
920  A 2D array of weight values that tapers smoothly to zero away from detected sources.
921  Set to a placeholder value of 1.0 if ``self.config.useModelWeights`` is False.
922  refImage : `lsst.afw.image.Image`
923  A reference image used to supply the default pixel values.
924  dcrWeights : `list` of `lsst.afw.image.Image`
925  Per-pixel weights for each subfilter.
926  Equal to 1/(number of unmasked images contributing to each pixel).
927 
928  Returns
929  -------
930  dcrModel : `lsst.pipe.tasks.DcrModel`
931  New model of the true sky after correcting chromatic effects.
932  """
933  newModelImages = []
934  for subfilter, model in enumerate(dcrModels):
935  residualsList = [next(residualGenerator) for residualGenerator in residualGeneratorList]
936  residual = np.sum(residualsList, axis=0)
937  residual *= dcrWeights[subfilter][dcrBBox].array
938  # `MaskedImage`s only support in-place addition, so rename for readability
939  newModel = model[dcrBBox].clone()
940  newModel.array += residual
941  # Catch any invalid values
942  badPixels = ~np.isfinite(newModel.array)
943  newModel.array[badPixels] = model[dcrBBox].array[badPixels]
944  if self.config.regularizeModelIterations > 0:
945  dcrModels.regularizeModelIter(subfilter, newModel, dcrBBox,
946  self.config.regularizeModelIterations,
947  self.config.regularizationWidth)
948  newModelImages.append(newModel)
949  if self.config.regularizeModelFrequency > 0:
950  dcrModels.regularizeModelFreq(newModelImages, dcrBBox, statsCtrl,
951  self.config.regularizeModelFrequency,
952  self.config.regularizationWidth)
953  dcrModels.conditionDcrModel(newModelImages, dcrBBox, gain=gain)
954  self.applyModelWeights(newModelImages, refImage[dcrBBox], modelWeights)
955  return DcrModel(newModelImages, dcrModels.filter, dcrModels.effectiveWavelength,
956  dcrModels.bandwidth, dcrModels.psf,
957  dcrModels.mask, dcrModels.variance)
958 
959  def calculateConvergence(self, dcrModels, subExposures, bbox, warpRefList, weightList, statsCtrl):
960  """Calculate a quality of fit metric for the matched templates.
961 
962  Parameters
963  ----------
964  dcrModels : `lsst.pipe.tasks.DcrModel`
965  Best fit model of the true sky after correcting chromatic effects.
966  subExposures : `dict` of `lsst.afw.image.ExposureF`
967  The pre-loaded exposures for the current subregion.
968  bbox : `lsst.geom.box.Box2I`
969  Sub-region to coadd
970  warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
971  `lsst.daf.persistence.ButlerDataRef`
972  The data references to the input warped exposures.
973  weightList : `list` of `float`
974  The weight to give each input exposure in the coadd
975  statsCtrl : `lsst.afw.math.StatisticsControl`
976  Statistics control object for coadd
977 
978  Returns
979  -------
980  convergenceMetric : `float`
981  Quality of fit metric for all input exposures, within the sub-region
982  """
983  significanceImage = np.abs(dcrModels.getReferenceImage(bbox))
984  nSigma = 3.
985  significanceImage += nSigma*dcrModels.calculateNoiseCutoff(dcrModels[1], statsCtrl,
986  bufferSize=self.bufferSize)
987  if np.max(significanceImage) == 0:
988  significanceImage += 1.
989  weight = 0
990  metric = 0.
991  metricList = {}
992  for warpExpRef, expWeight in zip(warpRefList, weightList):
993  visit = warpExpRef.dataId["visit"]
994  exposure = subExposures[visit][bbox]
995  singleMetric = self.calculateSingleConvergence(dcrModels, exposure, significanceImage, statsCtrl)
996  metric += singleMetric
997  metricList[visit] = singleMetric
998  weight += 1.
999  self.log.info("Individual metrics:\n%s", metricList)
1000  return 1.0 if weight == 0.0 else metric/weight
1001 
1002  def calculateSingleConvergence(self, dcrModels, exposure, significanceImage, statsCtrl):
1003  """Calculate a quality of fit metric for a single matched template.
1004 
1005  Parameters
1006  ----------
1007  dcrModels : `lsst.pipe.tasks.DcrModel`
1008  Best fit model of the true sky after correcting chromatic effects.
1009  exposure : `lsst.afw.image.ExposureF`
1010  The input warped exposure to evaluate.
1011  significanceImage : `numpy.ndarray`
1012  Array of weights for each pixel corresponding to its significance
1013  for the convergence calculation.
1014  statsCtrl : `lsst.afw.math.StatisticsControl`
1015  Statistics control object for coadd
1016 
1017  Returns
1018  -------
1019  convergenceMetric : `float`
1020  Quality of fit metric for one exposure, within the sub-region.
1021  """
1022  convergeMask = exposure.mask.getPlaneBitMask(self.config.convergenceMaskPlanes)
1023  templateImage = dcrModels.buildMatchedTemplate(exposure=exposure,
1024  order=self.config.imageInterpOrder,
1025  splitSubfilters=self.config.splitSubfilters,
1026  splitThreshold=self.config.splitThreshold,
1027  amplifyModel=self.config.accelerateModel)
1028  diffVals = np.abs(exposure.image.array - templateImage.array)*significanceImage
1029  refVals = np.abs(exposure.image.array + templateImage.array)*significanceImage/2.
1030 
1031  finitePixels = np.isfinite(diffVals)
1032  goodMaskPixels = (exposure.mask.array & statsCtrl.getAndMask()) == 0
1033  convergeMaskPixels = exposure.mask.array & convergeMask > 0
1034  usePixels = finitePixels & goodMaskPixels & convergeMaskPixels
1035  if np.sum(usePixels) == 0:
1036  metric = 0.
1037  else:
1038  diffUse = diffVals[usePixels]
1039  refUse = refVals[usePixels]
1040  metric = np.sum(diffUse/np.median(diffUse))/np.sum(refUse/np.median(diffUse))
1041  return metric
1042 
1043  def stackCoadd(self, dcrCoadds):
1044  """Add a list of sub-band coadds together.
1045 
1046  Parameters
1047  ----------
1048  dcrCoadds : `list` of `lsst.afw.image.ExposureF`
1049  A list of coadd exposures, each exposure containing
1050  the model for one subfilter.
1051 
1052  Returns
1053  -------
1054  coaddExposure : `lsst.afw.image.ExposureF`
1055  A single coadd exposure that is the sum of the sub-bands.
1056  """
1057  coaddExposure = dcrCoadds[0].clone()
1058  for coadd in dcrCoadds[1:]:
1059  coaddExposure.maskedImage += coadd.maskedImage
1060  return coaddExposure
1061 
1062  def fillCoadd(self, dcrModels, skyInfo, warpRefList, weightList, calibration=None, coaddInputs=None,
1063  mask=None, variance=None):
1064  """Create a list of coadd exposures from a list of masked images.
1065 
1066  Parameters
1067  ----------
1068  dcrModels : `lsst.pipe.tasks.DcrModel`
1069  Best fit model of the true sky after correcting chromatic effects.
1070  skyInfo : `lsst.pipe.base.Struct`
1071  Patch geometry information, from getSkyInfo
1072  warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
1073  `lsst.daf.persistence.ButlerDataRef`
1074  The data references to the input warped exposures.
1075  weightList : `list` of `float`
1076  The weight to give each input exposure in the coadd
1077  calibration : `lsst.afw.Image.PhotoCalib`, optional
1078  Scale factor to set the photometric calibration of an exposure.
1079  coaddInputs : `lsst.afw.Image.CoaddInputs`, optional
1080  A record of the observations that are included in the coadd.
1081  mask : `lsst.afw.image.Mask`, optional
1082  Optional mask to override the values in the final coadd.
1083  variance : `lsst.afw.image.Image`, optional
1084  Optional variance plane to override the values in the final coadd.
1085 
1086  Returns
1087  -------
1088  dcrCoadds : `list` of `lsst.afw.image.ExposureF`
1089  A list of coadd exposures, each exposure containing
1090  the model for one subfilter.
1091  """
1092  dcrCoadds = []
1093  refModel = dcrModels.getReferenceImage()
1094  for model in dcrModels:
1095  if self.config.accelerateModel > 1:
1096  model.array = (model.array - refModel)*self.config.accelerateModel + refModel
1097  coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs)
1098  if calibration is not None:
1099  coaddExposure.setPhotoCalib(calibration)
1100  if coaddInputs is not None:
1101  coaddExposure.getInfo().setCoaddInputs(coaddInputs)
1102  # Set the metadata for the coadd, including PSF and aperture corrections.
1103  self.assembleMetadata(coaddExposure, warpRefList, weightList)
1104  # Overwrite the PSF
1105  coaddExposure.setPsf(dcrModels.psf)
1106  coaddUtils.setCoaddEdgeBits(dcrModels.mask[skyInfo.bbox], dcrModels.variance[skyInfo.bbox])
1107  maskedImage = afwImage.MaskedImageF(dcrModels.bbox)
1108  maskedImage.image = model
1109  maskedImage.mask = dcrModels.mask
1110  maskedImage.variance = dcrModels.variance
1111  coaddExposure.setMaskedImage(maskedImage[skyInfo.bbox])
1112  coaddExposure.setPhotoCalib(self.scaleZeroPoint.getPhotoCalib())
1113  if mask is not None:
1114  coaddExposure.setMask(mask)
1115  if variance is not None:
1116  coaddExposure.setVariance(variance)
1117  dcrCoadds.append(coaddExposure)
1118  return dcrCoadds
1119 
1120  def calculateGain(self, convergenceList, gainList):
1121  """Calculate the gain to use for the current iteration.
1122 
1123  After calculating a new DcrModel, each value is averaged with the
1124  value in the corresponding pixel from the previous iteration. This
1125  reduces oscillating solutions that iterative techniques are plagued by,
1126  and speeds convergence. By far the biggest changes to the model
1127  happen in the first couple iterations, so we can also use a more
1128  aggressive gain later when the model is changing slowly.
1129 
1130  Parameters
1131  ----------
1132  convergenceList : `list` of `float`
1133  The quality of fit metric from each previous iteration.
1134  gainList : `list` of `float`
1135  The gains used in each previous iteration: appended with the new
1136  gain value.
1137  Gains are numbers between ``self.config.baseGain`` and 1.
1138 
1139  Returns
1140  -------
1141  gain : `float`
1142  Relative weight to give the new solution when updating the model.
1143  A value of 1.0 gives equal weight to both solutions.
1144 
1145  Raises
1146  ------
1147  ValueError
1148  If ``len(convergenceList) != len(gainList)+1``.
1149  """
1150  nIter = len(convergenceList)
1151  if nIter != len(gainList) + 1:
1152  raise ValueError("convergenceList (%d) must be one element longer than gainList (%d)."
1153  % (len(convergenceList), len(gainList)))
1154 
1155  if self.config.baseGain is None:
1156  # If ``baseGain`` is not set, calculate it from the number of DCR subfilters
1157  # The more subfilters being modeled, the lower the gain should be.
1158  baseGain = 1./(self.config.dcrNumSubfilters - 1)
1159  else:
1160  baseGain = self.config.baseGain
1161 
1162  if self.config.useProgressiveGain and nIter > 2:
1163  # To calculate the best gain to use, compare the past gains that have been used
1164  # with the resulting convergences to estimate the best gain to use.
1165  # Algorithmically, this is a Kalman filter.
1166  # If forward modeling proceeds perfectly, the convergence metric should
1167  # asymptotically approach a final value.
1168  # We can estimate that value from the measured changes in convergence
1169  # weighted by the gains used in each previous iteration.
1170  estFinalConv = [((1 + gainList[i])*convergenceList[i + 1] - convergenceList[i])/gainList[i]
1171  for i in range(nIter - 1)]
1172  # The convergence metric is strictly positive, so if the estimated final convergence is
1173  # less than zero, force it to zero.
1174  estFinalConv = np.array(estFinalConv)
1175  estFinalConv[estFinalConv < 0] = 0
1176  # Because the estimate may slowly change over time, only use the most recent measurements.
1177  estFinalConv = np.median(estFinalConv[max(nIter - 5, 0):])
1178  lastGain = gainList[-1]
1179  lastConv = convergenceList[-2]
1180  newConv = convergenceList[-1]
1181  # The predicted convergence is the value we would get if the new model calculated
1182  # in the previous iteration was perfect. Recall that the updated model that is
1183  # actually used is the gain-weighted average of the new and old model,
1184  # so the convergence would be similarly weighted.
1185  predictedConv = (estFinalConv*lastGain + lastConv)/(1. + lastGain)
1186  # If the measured and predicted convergence are very close, that indicates
1187  # that our forward model is accurate and we can use a more aggressive gain
1188  # If the measured convergence is significantly worse (or better!) than predicted,
1189  # that indicates that the model is not converging as expected and
1190  # we should use a more conservative gain.
1191  delta = (predictedConv - newConv)/((lastConv - estFinalConv)/(1 + lastGain))
1192  newGain = 1 - abs(delta)
1193  # Average the gains to prevent oscillating solutions.
1194  newGain = (newGain + lastGain)/2.
1195  gain = max(baseGain, newGain)
1196  else:
1197  gain = baseGain
1198  gainList.append(gain)
1199  return gain
1200 
1201  def calculateModelWeights(self, dcrModels, dcrBBox):
1202  """Build an array that smoothly tapers to 0 away from detected sources.
1203 
1204  Parameters
1205  ----------
1206  dcrModels : `lsst.pipe.tasks.DcrModel`
1207  Best fit model of the true sky after correcting chromatic effects.
1208  dcrBBox : `lsst.geom.box.Box2I`
1209  Sub-region of the coadd which includes a buffer to allow for DCR.
1210 
1211  Returns
1212  -------
1213  weights : `numpy.ndarray` or `float`
1214  A 2D array of weight values that tapers smoothly to zero away from detected sources.
1215  Set to a placeholder value of 1.0 if ``self.config.useModelWeights`` is False.
1216 
1217  Raises
1218  ------
1219  ValueError
1220  If ``useModelWeights`` is set and ``modelWeightsWidth`` is negative.
1221  """
1222  if not self.config.useModelWeights:
1223  return 1.0
1224  if self.config.modelWeightsWidth < 0:
1225  raise ValueError("modelWeightsWidth must not be negative if useModelWeights is set")
1226  convergeMask = dcrModels.mask.getPlaneBitMask(self.config.convergenceMaskPlanes)
1227  convergeMaskPixels = dcrModels.mask[dcrBBox].array & convergeMask > 0
1228  weights = np.zeros_like(dcrModels[0][dcrBBox].array)
1229  weights[convergeMaskPixels] = 1.
1230  weights = ndimage.filters.gaussian_filter(weights, self.config.modelWeightsWidth)
1231  weights /= np.max(weights)
1232  return weights
1233 
1234  def applyModelWeights(self, modelImages, refImage, modelWeights):
1235  """Smoothly replace model pixel values with those from a
1236  reference at locations away from detected sources.
1237 
1238  Parameters
1239  ----------
1240  modelImages : `list` of `lsst.afw.image.Image`
1241  The new DCR model images from the current iteration.
1242  The values will be modified in place.
1243  refImage : `lsst.afw.image.MaskedImage`
1244  A reference image used to supply the default pixel values.
1245  modelWeights : `numpy.ndarray` or `float`
1246  A 2D array of weight values that tapers smoothly to zero away from detected sources.
1247  Set to a placeholder value of 1.0 if ``self.config.useModelWeights`` is False.
1248  """
1249  if self.config.useModelWeights:
1250  for model in modelImages:
1251  model.array *= modelWeights
1252  model.array += refImage.array*(1. - modelWeights)/self.config.dcrNumSubfilters
1253 
1254  def loadSubExposures(self, bbox, statsCtrl, warpRefList, imageScalerList, spanSetMaskList):
1255  """Pre-load sub-regions of a list of exposures.
1256 
1257  Parameters
1258  ----------
1259  bbox : `lsst.geom.box.Box2I`
1260  Sub-region to coadd
1261  statsCtrl : `lsst.afw.math.StatisticsControl`
1262  Statistics control object for coadd
1263  warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
1264  `lsst.daf.persistence.ButlerDataRef`
1265  The data references to the input warped exposures.
1266  imageScalerList : `list` of `lsst.pipe.task.ImageScaler`
1267  The image scalars correct for the zero point of the exposures.
1268  spanSetMaskList : `list` of `dict` containing spanSet lists, or None
1269  Each element is dict with keys = mask plane name to add the spans to
1270 
1271  Returns
1272  -------
1273  subExposures : `dict`
1274  The `dict` keys are the visit IDs,
1275  and the values are `lsst.afw.image.ExposureF`
1276  The pre-loaded exposures for the current subregion.
1277  The variance plane contains weights, and not the variance
1278  """
1279  tempExpName = self.getTempExpDatasetName(self.warpType)
1280  zipIterables = zip(warpRefList, imageScalerList, spanSetMaskList)
1281  subExposures = {}
1282  for warpExpRef, imageScaler, altMaskSpans in zipIterables:
1283  if isinstance(warpExpRef, DeferredDatasetHandle):
1284  exposure = warpExpRef.get(parameters={'bbox': bbox})
1285  else:
1286  exposure = warpExpRef.get(tempExpName + "_sub", bbox=bbox)
1287  visit = warpExpRef.dataId["visit"]
1288  if altMaskSpans is not None:
1289  self.applyAltMaskPlanes(exposure.mask, altMaskSpans)
1290  imageScaler.scaleMaskedImage(exposure.maskedImage)
1291  # Note that the variance plane here is used to store weights, not the actual variance
1292  exposure.variance.array[:, :] = 0.
1293  # Set the weight of unmasked pixels to 1.
1294  exposure.variance.array[(exposure.mask.array & statsCtrl.getAndMask()) == 0] = 1.
1295  # Set the image value of masked pixels to zero.
1296  # This eliminates needing the mask plane when stacking images in ``newModelFromResidual``
1297  exposure.image.array[(exposure.mask.array & statsCtrl.getAndMask()) > 0] = 0.
1298  subExposures[visit] = exposure
1299  return subExposures
1300 
1301  def selectCoaddPsf(self, templateCoadd, warpRefList):
1302  """Compute the PSF of the coadd from the exposures with the best seeing.
1303 
1304  Parameters
1305  ----------
1306  templateCoadd : `lsst.afw.image.ExposureF`
1307  The initial coadd exposure before accounting for DCR.
1308  warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` or
1309  `lsst.daf.persistence.ButlerDataRef`
1310  The data references to the input warped exposures.
1311 
1312  Returns
1313  -------
1314  psf : `lsst.meas.algorithms.CoaddPsf`
1315  The average PSF of the input exposures with the best seeing.
1316  """
1317  sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
1318  tempExpName = self.getTempExpDatasetName(self.warpType)
1319  # Note: ``ccds`` is a `lsst.afw.table.ExposureCatalog` with one entry per ccd and per visit
1320  # If there are multiple ccds, it will have that many times more elements than ``warpExpRef``
1321  ccds = templateCoadd.getInfo().getCoaddInputs().ccds
1322  psfRefSize = templateCoadd.getPsf().computeShape().getDeterminantRadius()*sigma2fwhm
1323  psfSizes = np.zeros(len(ccds))
1324  ccdVisits = np.array(ccds["visit"])
1325  for warpExpRef in warpRefList:
1326  if isinstance(warpExpRef, DeferredDatasetHandle):
1327  # Gen 3 API
1328  psf = warpExpRef.get(component="psf")
1329  else:
1330  # Gen 2 API. Delete this when Gen 2 retired
1331  psf = warpExpRef.get(tempExpName).getPsf()
1332  visit = warpExpRef.dataId["visit"]
1333  psfSize = psf.computeShape().getDeterminantRadius()*sigma2fwhm
1334  psfSizes[ccdVisits == visit] = psfSize
1335  # Note that the input PSFs include DCR, which should be absent from the DcrCoadd
1336  # The selected PSFs are those that have a FWHM less than or equal to the smaller
1337  # of the mean or median FWHM of the input exposures.
1338  sizeThreshold = min(np.median(psfSizes), psfRefSize)
1339  goodPsfs = psfSizes <= sizeThreshold
1340  psf = measAlg.CoaddPsf(ccds[goodPsfs], templateCoadd.getWcs(),
1341  self.config.coaddPsf.makeControl())
1342  return psf
int min
int max
An integer coordinate rectangle.
Definition: Box.h:55
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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
Extent< int, N > ceil(Extent< double, N > const &input) noexcept
Return the component-wise ceil (round towards more positive).
Definition: Extent.cc:118
def applyDcr(image, dcr, useInverse=False, splitSubfilters=False, splitThreshold=0., doPrefilter=True, order=3)
Definition: dcrModel.py:694
def calculateDcr(visitInfo, wcs, effectiveWavelength, bandwidth, dcrNumSubfilters, splitSubfilters=False)
Definition: dcrModel.py:759
def run(self, coaddExposures, bbox, wcs)
Definition: getTemplate.py:603
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs)
def makeSkyInfo(skyMap, tractId, patchId)
Definition: coaddBase.py:289
def loadSubExposures(self, bbox, statsCtrl, warpRefList, imageScalerList, spanSetMaskList)
def fillCoadd(self, dcrModels, skyInfo, warpRefList, weightList, calibration=None, coaddInputs=None, mask=None, variance=None)
def applyModelWeights(self, modelImages, refImage, modelWeights)
def calculateSingleConvergence(self, dcrModels, exposure, significanceImage, statsCtrl)
def calculateConvergence(self, dcrModels, subExposures, bbox, warpRefList, weightList, statsCtrl)
def dcrAssembleSubregion(self, dcrModels, subExposures, bbox, dcrBBox, warpRefList, statsCtrl, convergenceMetric, gain, modelWeights, refImage, dcrWeights)
def calculateGain(self, convergenceList, gainList)
def calculateModelWeights(self, dcrModels, dcrBBox)
def newModelFromResidual(self, dcrModels, residualGeneratorList, dcrBBox, statsCtrl, gain, modelWeights, refImage, dcrWeights)
def selectCoaddPsf(self, templateCoadd, warpRefList)
def dcrResiduals(self, residual, visitInfo, wcs, effectiveWavelength, bandwidth)
Angle abs(Angle const &a)
Definition: Angle.h:106