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