LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
snapCombine.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010, 2011 LSST Corporation.
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 lsst.coadd.utils import addToCoadd, setCoaddEdgeBits
29 from lsst.ip.diffim import SnapPsfMatchTask
30 from lsst.meas.algorithms import SourceDetectionTask
31 from lsst.meas.base import SingleFrameMeasurementTask
32 import lsst.meas.algorithms as measAlg
33 
34 from .repair import RepairTask
35 from .calibrate import InitialPsfConfig
36 
37 class SnapCombineConfig(pexConfig.Config):
38  doRepair = pexConfig.Field(
39  dtype = bool,
40  doc = "Repair images (CR reject and interpolate) before combining",
41  default = True,
42  )
43  repairPsfFwhm = pexConfig.Field(
44  dtype = float,
45  doc = "Psf FWHM (pixels) used to detect CRs",
46  default = 2.5,
47  )
48  doDiffIm = pexConfig.Field(
49  dtype = bool,
50  doc = "Perform difference imaging before combining",
51  default = False,
52  )
53  doPsfMatch = pexConfig.Field(
54  dtype = bool,
55  doc = "Perform PSF matching for difference imaging (ignored if doDiffIm false)",
56  default = True,
57  )
58  doMeasurement = pexConfig.Field(
59  dtype = bool,
60  doc = "Measure difference sources (ignored if doDiffIm false)",
61  default = True,
62  )
63  badMaskPlanes = pexConfig.ListField(
64  dtype = str,
65  doc = "Mask planes that, if set, the associated pixels are not included in the combined exposure; "
66  "DETECTED excludes cosmic rays",
67  default = ("DETECTED",),
68  )
69  averageKeys = pexConfig.ListField(
70  dtype = str,
71  doc = "List of float metadata keys to average when combining snaps, e.g. float positions and dates; "
72  "non-float data must be handled by overriding the fixMetadata method",
73  optional = True,
74 
75  )
76  sumKeys = pexConfig.ListField(
77  dtype = str,
78  doc = "List of float or int metadata keys to sum when combining snaps, e.g. exposure time; "
79  "non-float, non-int data must be handled by overriding the fixMetadata method",
80  optional = True,
81  )
82 
83  repair = pexConfig.ConfigurableField(target = RepairTask, doc = "")
84  diffim = pexConfig.ConfigurableField(target = SnapPsfMatchTask, doc = "")
85  detection = pexConfig.ConfigurableField(target = SourceDetectionTask, doc = "")
86  initialPsf = pexConfig.ConfigField(dtype = InitialPsfConfig, doc = "")
87  measurement = pexConfig.ConfigurableField(target = SingleFrameMeasurementTask, doc = "")
88 
89  def setDefaults(self):
90  self.detection.thresholdPolarity = "both"
91 
92  def validate(self):
93  if self.detection.thresholdPolarity != "both":
94  raise ValueError("detection.thresholdPolarity must be 'both' for SnapCombineTask")
95 
96 class SnapCombineTask(pipeBase.Task):
97  ConfigClass = SnapCombineConfig
98  _DefaultName = "snapCombine"
99 
100  def __init__(self, *args, **kwargs):
101  pipeBase.Task.__init__(self, *args, **kwargs)
102  self.makeSubtask("repair")
103  self.makeSubtask("diffim")
104  self.schema = afwTable.SourceTable.makeMinimalSchema()
106  self.makeSubtask("detection", schema=self.schema)
107  if self.config.doMeasurement:
108  self.makeSubtask("measurement", schema=self.schema, algMetadata=self.algMetadata)
109 
110  @pipeBase.timeMethod
111  def run(self, snap0, snap1, defects=None):
112  """Combine two snaps
113 
114  @param[in] snap0: snapshot exposure 0
115  @param[in] snap1: snapshot exposure 1
116  @defects[in] defect list (for repair task)
117  @return a pipe_base Struct with fields:
118  - exposure: snap-combined exposure
119  - sources: detected sources, or None if detection not performed
120  """
121  # initialize optional outputs
122  sources = None
123 
124  if self.config.doRepair:
125  self.log.info("snapCombine repair")
126  psf = self.makeInitialPsf(snap0, fwhmPix=self.config.repairPsfFwhm)
127  snap0.setPsf(psf)
128  snap1.setPsf(psf)
129  self.repair.run(snap0, defects=defects, keepCRs=False)
130  self.repair.run(snap1, defects=defects, keepCRs=False)
131  self.display('repair0', exposure=snap0)
132  self.display('repair1', exposure=snap1)
133 
134  if self.config.doDiffIm:
135  if self.config.doPsfMatch:
136  self.log.info("snapCombine psfMatch")
137  diffRet = self.diffim.run(snap0, snap1, "subtractExposures")
138  diffExp = diffRet.subtractedImage
139 
140  # Measure centroid and width of kernel; dependent on ticket #1980
141  # Useful diagnostic for the degree of astrometric shift between snaps.
142  diffKern = diffRet.psfMatchingKernel
143  width, height = diffKern.getDimensions()
144  # TBD...
145  #psfAttr = measAlg.PsfAttributes(diffKern, width//2, height//2)
146 
147  else:
148  diffExp = afwImage.ExposureF(snap0, True)
149  diffMi = diffExp.getMaskedImage()
150  diffMi -= snap1.getMaskedImage()
151 
152  psf = self.makeInitialPsf(snap0)
153  diffExp.setPsf(psf)
154  table = afwTable.SourceTable.make(self.schema)
155  table.setMetadata(self.algMetadata)
156  detRet = self.detection.makeSourceCatalog(table, diffExp)
157  sources = detRet.sources
158  fpSets = detRet.fpSets
159  if self.config.doMeasurement:
160  self.measurement.measure(diffExp, sources)
161 
162  mask0 = snap0.getMaskedImage().getMask()
163  mask1 = snap1.getMaskedImage().getMask()
164  fpSets.positive.setMask(mask0, "DETECTED")
165  fpSets.negative.setMask(mask1, "DETECTED")
166 
167  maskD = diffExp.getMaskedImage().getMask()
168  fpSets.positive.setMask(maskD, "DETECTED")
169  fpSets.negative.setMask(maskD, "DETECTED_NEGATIVE")
170 
171  combinedExp = self.addSnaps(snap0, snap1)
172 
173  return pipeBase.Struct(
174  exposure = combinedExp,
175  sources = sources,
176  )
177 
178  def addSnaps(self, snap0, snap1):
179  """Add two snap exposures together, returning a new exposure
180 
181  @param[in] snap0 snap exposure 0
182  @param[in] snap1 snap exposure 1
183  @return combined exposure
184  """
185  self.log.info("snapCombine addSnaps")
186 
187  combinedExp = snap0.Factory(snap0, True)
188  combinedMi = combinedExp.getMaskedImage()
189  combinedMi.set(0)
190 
191  weightMap = combinedMi.getImage().Factory(combinedMi.getBBox())
192  weight = 1.0
193  badPixelMask = afwImage.MaskU.getPlaneBitMask(self.config.badMaskPlanes)
194  addToCoadd(combinedMi, weightMap, snap0.getMaskedImage(), badPixelMask, weight)
195  addToCoadd(combinedMi, weightMap, snap1.getMaskedImage(), badPixelMask, weight)
196 
197  # pre-scaling the weight map instead of post-scaling the combinedMi saves a bit of time
198  # because the weight map is a simple Image instead of a MaskedImage
199  weightMap *= 0.5 # so result is sum of both images, instead of average
200  combinedMi /= weightMap
201  setCoaddEdgeBits(combinedMi.getMask(), weightMap)
202 
203  # note: none of the inputs has a valid Calib object, so that is not touched
204  # Filter was already copied
205 
206  combinedMetadata = combinedExp.getMetadata()
207  metadata0 = snap0.getMetadata()
208  metadata1 = snap1.getMetadata()
209  self.fixMetadata(combinedMetadata, metadata0, metadata1)
210 
211  return combinedExp
212 
213  def fixMetadata(self, combinedMetadata, metadata0, metadata1):
214  """Fix the metadata of the combined exposure (in place)
215 
216  This implementation handles items specified by config.averageKeys and config.sumKeys,
217  which have data type restrictions. To handle other data types (such as sexagesimal
218  positions and ISO dates) you must supplement this method with your own code.
219 
220  @param[in,out] combinedMetadata metadata of combined exposure;
221  on input this is a deep copy of metadata0 (a PropertySet)
222  @param[in] metadata0 metadata of snap0 (a PropertySet)
223  @param[in] metadata1 metadata of snap1 (a PropertySet)
224 
225  @note the inputs are presently PropertySets due to ticket #2542. However, in some sense
226  they are just PropertyLists that are missing some methods. In particular: comments and order
227  are preserved if you alter an existing value with set(key, value).
228  """
229  keyDoAvgList = []
230  if self.config.averageKeys:
231  keyDoAvgList += [(key, 1) for key in self.config.averageKeys]
232  if self.config.sumKeys:
233  keyDoAvgList += [(key, 0) for key in self.config.sumKeys]
234  for key, doAvg in keyDoAvgList:
235  opStr = "average" if doAvg else "sum"
236  try:
237  val0 = metadata0.get(key)
238  val1 = metadata1.get(key)
239  except Exception:
240  self.log.warn("Could not %s metadata %r: missing from one or both exposures" % (opStr, key,))
241  continue
242 
243  try:
244  combinedVal = val0 + val1
245  if doAvg:
246  combinedVal /= 2.0
247  except Exception:
248  self.log.warn("Could not %s metadata %r: value %r and/or %r not numeric" % \
249  (opStr, key, val0, val1))
250  continue
251 
252  combinedMetadata.set(key, combinedVal)
253 
254  def makeInitialPsf(self, exposure, fwhmPix=None):
255  """Initialise the detection procedure by setting the PSF and WCS
256 
257  @param exposure Exposure to process
258  @return PSF, WCS
259  """
260  assert exposure, "No exposure provided"
261  wcs = exposure.getWcs()
262  assert wcs, "No wcs in exposure"
263 
264  if fwhmPix is None:
265  fwhmPix = self.config.initialPsf.fwhm / wcs.pixelScale().asArcseconds()
266 
267  size = self.config.initialPsf.size
268  model = self.config.initialPsf.model
269  self.log.info("installInitialPsf fwhm=%s pixels; size=%s pixels" % (fwhmPix, size))
270  psfCls = getattr(measAlg, model + "Psf")
271  psf = psfCls(size, size, fwhmPix/(2.0*num.sqrt(2*num.log(2.0))))
272  return psf
Class for storing ordered metadata with comments.
Definition: PropertyList.h:81
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
lsst::afw::geom::Box2I addToCoadd(lsst::afw::image::Image< CoaddPixelT > &coadd, lsst::afw::image::Image< WeightPixelT > &weightMap, lsst::afw::image::Image< CoaddPixelT > const &image, WeightPixelT weight)
add good pixels from an image to a coadd and associated weight map
Definition: addToCoadd.cc:126