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
snapCombine.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2016 AURA/LSST.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
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 <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 import numpy as num
23 import lsst.pex.config as pexConfig
24 import lsst.daf.base as dafBase
25 import lsst.afw.image as afwImage
26 import lsst.afw.table as afwTable
27 import lsst.pipe.base as pipeBase
28 from lsstDebug import getDebugFrame
29 from lsst.afw.display import getDisplay
30 from lsst.coadd.utils import addToCoadd, setCoaddEdgeBits
31 from lsst.ip.diffim import SnapPsfMatchTask
32 from lsst.meas.algorithms import SourceDetectionTask
33 from lsst.meas.base import SingleFrameMeasurementTask
34 import lsst.meas.algorithms as measAlg
35 from lsst.utils.timer import timeMethod
36 
37 from .repair import RepairTask
38 
39 
40 class InitialPsfConfig(pexConfig.Config):
41  """!Describes the initial PSF used for detection and measurement before we do PSF determination."""
42 
43  model = pexConfig.ChoiceField(
44  dtype=str,
45  doc="PSF model type",
46  default="SingleGaussian",
47  allowed={
48  "SingleGaussian": "Single Gaussian model",
49  "DoubleGaussian": "Double Gaussian model",
50  },
51  )
52  pixelScale = pexConfig.Field(
53  dtype=float,
54  doc="Pixel size (arcsec). Only needed if no Wcs is provided",
55  default=0.25,
56  )
57  fwhm = pexConfig.Field(
58  dtype=float,
59  doc="FWHM of PSF model (arcsec)",
60  default=1.0,
61  )
62  size = pexConfig.Field(
63  dtype=int,
64  doc="Size of PSF model (pixels)",
65  default=15,
66  )
67 
68 
69 class SnapCombineConfig(pexConfig.Config):
70  doRepair = pexConfig.Field(
71  dtype=bool,
72  doc="Repair images (CR reject and interpolate) before combining",
73  default=True,
74  )
75  repairPsfFwhm = pexConfig.Field(
76  dtype=float,
77  doc="Psf FWHM (pixels) used to detect CRs",
78  default=2.5,
79  )
80  doDiffIm = pexConfig.Field(
81  dtype=bool,
82  doc="Perform difference imaging before combining",
83  default=False,
84  )
85  doPsfMatch = pexConfig.Field(
86  dtype=bool,
87  doc="Perform PSF matching for difference imaging (ignored if doDiffIm false)",
88  default=True,
89  )
90  doMeasurement = pexConfig.Field(
91  dtype=bool,
92  doc="Measure difference sources (ignored if doDiffIm false)",
93  default=True,
94  )
95  badMaskPlanes = pexConfig.ListField(
96  dtype=str,
97  doc="Mask planes that, if set, the associated pixels are not included in the combined exposure; "
98  "DETECTED excludes cosmic rays",
99  default=("DETECTED",),
100  )
101  averageKeys = pexConfig.ListField(
102  dtype=str,
103  doc="List of float metadata keys to average when combining snaps, e.g. float positions and dates; "
104  "non-float data must be handled by overriding the fixMetadata method",
105  optional=True,
106 
107  )
108  sumKeys = pexConfig.ListField(
109  dtype=str,
110  doc="List of float or int metadata keys to sum when combining snaps, e.g. exposure time; "
111  "non-float, non-int data must be handled by overriding the fixMetadata method",
112  optional=True,
113  )
114 
115  repair = pexConfig.ConfigurableField(target=RepairTask, doc="")
116  diffim = pexConfig.ConfigurableField(target=SnapPsfMatchTask, doc="")
117  detection = pexConfig.ConfigurableField(target=SourceDetectionTask, doc="")
118  initialPsf = pexConfig.ConfigField(dtype=InitialPsfConfig, doc="")
119  measurement = pexConfig.ConfigurableField(target=SingleFrameMeasurementTask, doc="")
120 
121  def setDefaults(self):
122  self.detectiondetection.thresholdPolarity = "both"
123 
124  def validate(self):
125  if self.detectiondetection.thresholdPolarity != "both":
126  raise ValueError("detection.thresholdPolarity must be 'both' for SnapCombineTask")
127 
128 
134 
135 
136 class SnapCombineTask(pipeBase.Task):
137  r"""!
138  \anchor SnapCombineTask_
139 
140  \brief Combine snaps.
141 
142  \section pipe_tasks_snapcombine_Contents Contents
143 
144  - \ref pipe_tasks_snapcombine_Debug
145 
146  \section pipe_tasks_snapcombine_Debug Debug variables
147 
148  The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a
149  flag \c -d to import \b debug.py from your \c PYTHONPATH; see <a
150  href="https://developer.lsst.io/stack/debug.html">Debugging Tasks with lsstDebug</a> for more
151  about \b debug.py files.
152 
153  The available variables in SnapCombineTask are:
154  <DL>
155  <DT> \c display
156  <DD> A dictionary containing debug point names as keys with frame number as value. Valid keys are:
157  <DL>
158  <DT> repair0
159  <DD> Display the first snap after repairing.
160  <DT> repair1
161  <DD> Display the second snap after repairing.
162  </DL>
163  </DD>
164  </DL>
165  """
166  ConfigClass = SnapCombineConfig
167  _DefaultName = "snapCombine"
168 
169  def __init__(self, *args, **kwargs):
170  pipeBase.Task.__init__(self, *args, **kwargs)
171  self.makeSubtask("repair")
172  self.makeSubtask("diffim")
173  self.schemaschema = afwTable.SourceTable.makeMinimalSchema()
174  self.algMetadataalgMetadata = dafBase.PropertyList()
175  self.makeSubtask("detection", schema=self.schemaschema)
176  if self.config.doMeasurement:
177  self.makeSubtask("measurement", schema=self.schemaschema, algMetadata=self.algMetadataalgMetadata)
178 
179  @timeMethod
180  def run(self, snap0, snap1, defects=None):
181  """Combine two snaps
182 
183  @param[in] snap0: snapshot exposure 0
184  @param[in] snap1: snapshot exposure 1
185  @defects[in] defect list (for repair task)
186  @return a pipe_base Struct with fields:
187  - exposure: snap-combined exposure
188  - sources: detected sources, or None if detection not performed
189  """
190  # initialize optional outputs
191  sources = None
192 
193  if self.config.doRepair:
194  self.log.info("snapCombine repair")
195  psf = self.makeInitialPsfmakeInitialPsf(snap0, fwhmPix=self.config.repairPsfFwhm)
196  snap0.setPsf(psf)
197  snap1.setPsf(psf)
198  self.repair.run(snap0, defects=defects, keepCRs=False)
199  self.repair.run(snap1, defects=defects, keepCRs=False)
200 
201  repair0frame = getDebugFrame(self._display, "repair0")
202  if repair0frame:
203  getDisplay(repair0frame).mtv(snap0)
204  repair1frame = getDebugFrame(self._display, "repair1")
205  if repair1frame:
206  getDisplay(repair1frame).mtv(snap1)
207 
208  if self.config.doDiffIm:
209  if self.config.doPsfMatch:
210  self.log.info("snapCombine psfMatch")
211  diffRet = self.diffim.run(snap0, snap1, "subtractExposures")
212  diffExp = diffRet.subtractedImage
213 
214  # Measure centroid and width of kernel; dependent on ticket #1980
215  # Useful diagnostic for the degree of astrometric shift between snaps.
216  diffKern = diffRet.psfMatchingKernel
217  width, height = diffKern.getDimensions()
218 
219  else:
220  diffExp = afwImage.ExposureF(snap0, True)
221  diffMi = diffExp.getMaskedImage()
222  diffMi -= snap1.getMaskedImage()
223 
224  psf = self.makeInitialPsfmakeInitialPsf(snap0)
225  diffExp.setPsf(psf)
226  table = afwTable.SourceTable.make(self.schemaschema)
227  table.setMetadata(self.algMetadataalgMetadata)
228  detRet = self.detection.run(table, diffExp)
229  sources = detRet.sources
230  fpSets = detRet.fpSets
231  if self.config.doMeasurement:
232  self.measurement.measure(diffExp, sources)
233 
234  mask0 = snap0.getMaskedImage().getMask()
235  mask1 = snap1.getMaskedImage().getMask()
236  fpSets.positive.setMask(mask0, "DETECTED")
237  fpSets.negative.setMask(mask1, "DETECTED")
238 
239  maskD = diffExp.getMaskedImage().getMask()
240  fpSets.positive.setMask(maskD, "DETECTED")
241  fpSets.negative.setMask(maskD, "DETECTED_NEGATIVE")
242 
243  combinedExp = self.addSnapsaddSnaps(snap0, snap1)
244 
245  return pipeBase.Struct(
246  exposure=combinedExp,
247  sources=sources,
248  )
249 
250  def addSnaps(self, snap0, snap1):
251  """Add two snap exposures together, returning a new exposure
252 
253  @param[in] snap0 snap exposure 0
254  @param[in] snap1 snap exposure 1
255  @return combined exposure
256  """
257  self.log.info("snapCombine addSnaps")
258 
259  combinedExp = snap0.Factory(snap0, True)
260  combinedMi = combinedExp.getMaskedImage()
261  combinedMi.set(0)
262 
263  weightMap = combinedMi.getImage().Factory(combinedMi.getBBox())
264  weight = 1.0
265  badPixelMask = afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes)
266  addToCoadd(combinedMi, weightMap, snap0.getMaskedImage(), badPixelMask, weight)
267  addToCoadd(combinedMi, weightMap, snap1.getMaskedImage(), badPixelMask, weight)
268 
269  # pre-scaling the weight map instead of post-scaling the combinedMi saves a bit of time
270  # because the weight map is a simple Image instead of a MaskedImage
271  weightMap *= 0.5 # so result is sum of both images, instead of average
272  combinedMi /= weightMap
273  setCoaddEdgeBits(combinedMi.getMask(), weightMap)
274 
275  # note: none of the inputs has a valid PhotoCalib object, so that is not touched
276  # Filter was already copied
277 
278  combinedMetadata = combinedExp.getMetadata()
279  metadata0 = snap0.getMetadata()
280  metadata1 = snap1.getMetadata()
281  self.fixMetadatafixMetadata(combinedMetadata, metadata0, metadata1)
282 
283  return combinedExp
284 
285  def fixMetadata(self, combinedMetadata, metadata0, metadata1):
286  """Fix the metadata of the combined exposure (in place)
287 
288  This implementation handles items specified by config.averageKeys and config.sumKeys,
289  which have data type restrictions. To handle other data types (such as sexagesimal
290  positions and ISO dates) you must supplement this method with your own code.
291 
292  @param[in,out] combinedMetadata metadata of combined exposure;
293  on input this is a deep copy of metadata0 (a PropertySet)
294  @param[in] metadata0 metadata of snap0 (a PropertySet)
295  @param[in] metadata1 metadata of snap1 (a PropertySet)
296 
297  @note the inputs are presently PropertySets due to ticket #2542. However, in some sense
298  they are just PropertyLists that are missing some methods. In particular: comments and order
299  are preserved if you alter an existing value with set(key, value).
300  """
301  keyDoAvgList = []
302  if self.config.averageKeys:
303  keyDoAvgList += [(key, 1) for key in self.config.averageKeys]
304  if self.config.sumKeys:
305  keyDoAvgList += [(key, 0) for key in self.config.sumKeys]
306  for key, doAvg in keyDoAvgList:
307  opStr = "average" if doAvg else "sum"
308  try:
309  val0 = metadata0.getScalar(key)
310  val1 = metadata1.getScalar(key)
311  except Exception:
312  self.log.warning("Could not %s metadata %r: missing from one or both exposures", opStr, key)
313  continue
314 
315  try:
316  combinedVal = val0 + val1
317  if doAvg:
318  combinedVal /= 2.0
319  except Exception:
320  self.log.warning("Could not %s metadata %r: value %r and/or %r not numeric",
321  opStr, key, val0, val1)
322  continue
323 
324  combinedMetadata.set(key, combinedVal)
325 
326  def makeInitialPsf(self, exposure, fwhmPix=None):
327  """Initialise the detection procedure by setting the PSF and WCS
328 
329  @param exposure Exposure to process
330  @return PSF, WCS
331  """
332  assert exposure, "No exposure provided"
333  wcs = exposure.getWcs()
334  assert wcs, "No wcs in exposure"
335 
336  if fwhmPix is None:
337  fwhmPix = self.config.initialPsf.fwhm / wcs.getPixelScale().asArcseconds()
338 
339  size = self.config.initialPsf.size
340  model = self.config.initialPsf.model
341  self.log.info("installInitialPsf fwhm=%s pixels; size=%s pixels", fwhmPix, size)
342  psfCls = getattr(measAlg, model + "Psf")
343  psf = psfCls(size, size, fwhmPix/(2.0*num.sqrt(2*num.log(2.0))))
344  return psf
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
Describes the initial PSF used for detection and measurement before we do PSF determination.
Definition: snapCombine.py:40
def fixMetadata(self, combinedMetadata, metadata0, metadata1)
Definition: snapCombine.py:285
def makeInitialPsf(self, exposure, fwhmPix=None)
Definition: snapCombine.py:326
def run(self, snap0, snap1, defects=None)
Definition: snapCombine.py:180
def __init__(self, *args, **kwargs)
Definition: snapCombine.py:169
def mtv(data, frame=None, title="", wcs=None, *args, **kwargs)
Definition: ds9.py:92
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
def measure(mi, x, y, size, statistic, stats)
Definition: fringe.py:517
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
def getDebugFrame(debugDisplay, name)
Definition: lsstDebug.py:95