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
imageDifference.py
Go to the documentation of this file.
1 # This file is part of pipe_tasks.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <https://www.gnu.org/licenses/>.
21 
22 import math
23 import random
24 import numpy
25 
26 import lsst.utils
27 import lsst.pex.config as pexConfig
28 import lsst.pipe.base as pipeBase
29 import lsst.daf.base as dafBase
30 import lsst.geom as geom
31 import lsst.afw.math as afwMath
32 import lsst.afw.table as afwTable
33 import lsst.meas.extensions.trailedSources # noqa: F401
34 from lsst.meas.astrom import AstrometryConfig, AstrometryTask
35 from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask
36 from lsst.meas.algorithms import LoadIndexedReferenceObjectsTask, SkyObjectsTask
37 from lsst.pipe.tasks.registerImage import RegisterTask
38 from lsst.pipe.tasks.scaleVariance import ScaleVarianceTask
39 from lsst.meas.algorithms import SourceDetectionTask, SingleGaussianPsf, ObjectSizeStarSelectorTask
40 from lsst.ip.diffim import (DipoleAnalysis, SourceFlagChecker, KernelCandidateF, makeKernelBasisList,
41  KernelCandidateQa, DiaCatalogSourceSelectorTask, DiaCatalogSourceSelectorConfig,
42  GetCoaddAsTemplateTask, GetCalexpAsTemplateTask, DipoleFitTask,
43  DecorrelateALKernelSpatialTask, subtractAlgorithmRegistry)
44 import lsst.ip.diffim.diffimTools as diffimTools
45 import lsst.ip.diffim.utils as diUtils
46 import lsst.afw.display as afwDisplay
47 from lsst.skymap import BaseSkyMap
48 from lsst.obs.base import ExposureIdInfo
49 from lsst.utils.timer import timeMethod
50 
51 __all__ = ["ImageDifferenceConfig", "ImageDifferenceTask"]
52 FwhmPerSigma = 2*math.sqrt(2*math.log(2))
53 IqrToSigma = 0.741
54 
55 
56 class ImageDifferenceTaskConnections(pipeBase.PipelineTaskConnections,
57  dimensions=("instrument", "visit", "detector", "skymap"),
58  defaultTemplates={"coaddName": "deep",
59  "skyMapName": "deep",
60  "warpTypeSuffix": "",
61  "fakesType": ""}):
62 
63  exposure = pipeBase.connectionTypes.Input(
64  doc="Input science exposure to subtract from.",
65  dimensions=("instrument", "visit", "detector"),
66  storageClass="ExposureF",
67  name="{fakesType}calexp"
68  )
69 
70  # TODO DM-22953
71  # kernelSources = pipeBase.connectionTypes.Input(
72  # doc="Source catalog produced in calibrate task for kernel candidate sources",
73  # name="src",
74  # storageClass="SourceCatalog",
75  # dimensions=("instrument", "visit", "detector"),
76  # )
77 
78  skyMap = pipeBase.connectionTypes.Input(
79  doc="Input definition of geometry/bbox and projection/wcs for template exposures",
80  name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
81  dimensions=("skymap", ),
82  storageClass="SkyMap",
83  )
84  coaddExposures = pipeBase.connectionTypes.Input(
85  doc="Input template to match and subtract from the exposure",
86  dimensions=("tract", "patch", "skymap", "band"),
87  storageClass="ExposureF",
88  name="{fakesType}{coaddName}Coadd{warpTypeSuffix}",
89  multiple=True,
90  deferLoad=True
91  )
92  dcrCoadds = pipeBase.connectionTypes.Input(
93  doc="Input DCR template to match and subtract from the exposure",
94  name="{fakesType}dcrCoadd{warpTypeSuffix}",
95  storageClass="ExposureF",
96  dimensions=("tract", "patch", "skymap", "band", "subfilter"),
97  multiple=True,
98  deferLoad=True
99  )
100  outputSchema = pipeBase.connectionTypes.InitOutput(
101  doc="Schema (as an example catalog) for output DIASource catalog.",
102  storageClass="SourceCatalog",
103  name="{fakesType}{coaddName}Diff_diaSrc_schema",
104  )
105  subtractedExposure = pipeBase.connectionTypes.Output(
106  doc="Output AL difference or Zogy proper difference image",
107  dimensions=("instrument", "visit", "detector"),
108  storageClass="ExposureF",
109  name="{fakesType}{coaddName}Diff_differenceExp",
110  )
111  scoreExposure = pipeBase.connectionTypes.Output(
112  doc="Output AL likelihood or Zogy score image",
113  dimensions=("instrument", "visit", "detector"),
114  storageClass="ExposureF",
115  name="{fakesType}{coaddName}Diff_scoreExp",
116  )
117  warpedExposure = pipeBase.connectionTypes.Output(
118  doc="Warped template used to create `subtractedExposure`.",
119  dimensions=("instrument", "visit", "detector"),
120  storageClass="ExposureF",
121  name="{fakesType}{coaddName}Diff_warpedExp",
122  )
123  matchedExposure = pipeBase.connectionTypes.Output(
124  doc="Warped template used to create `subtractedExposure`.",
125  dimensions=("instrument", "visit", "detector"),
126  storageClass="ExposureF",
127  name="{fakesType}{coaddName}Diff_matchedExp",
128  )
129  diaSources = pipeBase.connectionTypes.Output(
130  doc="Output detected diaSources on the difference image",
131  dimensions=("instrument", "visit", "detector"),
132  storageClass="SourceCatalog",
133  name="{fakesType}{coaddName}Diff_diaSrc",
134  )
135 
136  def __init__(self, *, config=None):
137  super().__init__(config=config)
138  if config.coaddName == 'dcr':
139  self.inputs.remove("coaddExposures")
140  else:
141  self.inputs.remove("dcrCoadds")
142  if not config.doWriteSubtractedExp:
143  self.outputs.remove("subtractedExposure")
144  if not config.doWriteScoreExp:
145  self.outputs.remove("scoreExposure")
146  if not config.doWriteWarpedExp:
147  self.outputs.remove("warpedExposure")
148  if not config.doWriteMatchedExp:
149  self.outputs.remove("matchedExposure")
150  if not config.doWriteSources:
151  self.outputs.remove("diaSources")
152 
153  # TODO DM-22953: Add support for refObjLoader (kernelSourcesFromRef)
154  # Make kernelSources optional
155 
156 
157 class ImageDifferenceConfig(pipeBase.PipelineTaskConfig,
158  pipelineConnections=ImageDifferenceTaskConnections):
159  """Config for ImageDifferenceTask
160  """
161  doAddCalexpBackground = pexConfig.Field(dtype=bool, default=False,
162  doc="Add background to calexp before processing it. "
163  "Useful as ipDiffim does background matching.")
164  doUseRegister = pexConfig.Field(dtype=bool, default=False,
165  doc="Re-compute astrometry on the template. "
166  "Use image-to-image registration to align template with "
167  "science image (AL only).")
168  doDebugRegister = pexConfig.Field(dtype=bool, default=False,
169  doc="Writing debugging data for doUseRegister")
170  doSelectSources = pexConfig.Field(dtype=bool, default=False,
171  doc="Select stars to use for kernel fitting (AL only)")
172  doSelectDcrCatalog = pexConfig.Field(dtype=bool, default=False,
173  doc="Select stars of extreme color as part "
174  "of the control sample (AL only)")
175  doSelectVariableCatalog = pexConfig.Field(dtype=bool, default=False,
176  doc="Select stars that are variable to be part "
177  "of the control sample (AL only)")
178  doSubtract = pexConfig.Field(dtype=bool, default=True, doc="Compute subtracted exposure?")
179  doPreConvolve = pexConfig.Field(dtype=bool, default=False,
180  doc="Not in use. Superseded by useScoreImageDetection.",
181  deprecated="This option superseded by useScoreImageDetection."
182  " Will be removed after v22.")
183  useScoreImageDetection = pexConfig.Field(
184  dtype=bool, default=False, doc="Calculate the pre-convolved AL likelihood or "
185  "the Zogy score image. Use it for source detection (if doDetection=True).")
186  doWriteScoreExp = pexConfig.Field(
187  dtype=bool, default=False, doc="Write AL likelihood or Zogy score exposure?")
188  doScaleTemplateVariance = pexConfig.Field(dtype=bool, default=False,
189  doc="Scale variance of the template before PSF matching")
190  doScaleDiffimVariance = pexConfig.Field(dtype=bool, default=True,
191  doc="Scale variance of the diffim before PSF matching. "
192  "You may do either this or template variance scaling, "
193  "or neither. (Doing both is a waste of CPU.)")
194  useGaussianForPreConvolution = pexConfig.Field(
195  dtype=bool, default=False, doc="Use a simple gaussian PSF model for pre-convolution "
196  "(oherwise use exposure PSF)? (AL and if useScoreImageDetection=True only)")
197  doDetection = pexConfig.Field(dtype=bool, default=True, doc="Detect sources?")
198  doDecorrelation = pexConfig.Field(dtype=bool, default=True,
199  doc="Perform diffim decorrelation to undo pixel correlation due to A&L "
200  "kernel convolution (AL only)? If True, also update the diffim PSF.")
201  doMerge = pexConfig.Field(dtype=bool, default=True,
202  doc="Merge positive and negative diaSources with grow radius "
203  "set by growFootprint")
204  doMatchSources = pexConfig.Field(dtype=bool, default=False,
205  doc="Match diaSources with input calexp sources and ref catalog sources")
206  doMeasurement = pexConfig.Field(dtype=bool, default=True, doc="Measure diaSources?")
207  doDipoleFitting = pexConfig.Field(dtype=bool, default=True, doc="Measure dipoles using new algorithm?")
208  doForcedMeasurement = pexConfig.Field(
209  dtype=bool,
210  default=True,
211  doc="Force photometer diaSource locations on PVI?")
212  doWriteSubtractedExp = pexConfig.Field(
213  dtype=bool, default=True, doc="Write difference exposure (AL and Zogy) ?")
214  doWriteWarpedExp = pexConfig.Field(
215  dtype=bool, default=False, doc="Write WCS, warped template coadd exposure?")
216  doWriteMatchedExp = pexConfig.Field(dtype=bool, default=False,
217  doc="Write warped and PSF-matched template coadd exposure?")
218  doWriteSources = pexConfig.Field(dtype=bool, default=True, doc="Write sources?")
219  doAddMetrics = pexConfig.Field(dtype=bool, default=False,
220  doc="Add columns to the source table to hold analysis metrics?")
221 
222  coaddName = pexConfig.Field(
223  doc="coadd name: typically one of deep, goodSeeing, or dcr",
224  dtype=str,
225  default="deep",
226  )
227  convolveTemplate = pexConfig.Field(
228  doc="Which image gets convolved (default = template)",
229  dtype=bool,
230  default=True
231  )
232  refObjLoader = pexConfig.ConfigurableField(
233  target=LoadIndexedReferenceObjectsTask,
234  doc="reference object loader",
235  )
236  astrometer = pexConfig.ConfigurableField(
237  target=AstrometryTask,
238  doc="astrometry task; used to match sources to reference objects, but not to fit a WCS",
239  )
240  sourceSelector = pexConfig.ConfigurableField(
241  target=ObjectSizeStarSelectorTask,
242  doc="Source selection algorithm",
243  )
244  subtract = subtractAlgorithmRegistry.makeField("Subtraction Algorithm", default="al")
245  decorrelate = pexConfig.ConfigurableField(
246  target=DecorrelateALKernelSpatialTask,
247  doc="Decorrelate effects of A&L kernel convolution on image difference, only if doSubtract is True. "
248  "If this option is enabled, then detection.thresholdValue should be set to 5.0 (rather than the "
249  "default of 5.5).",
250  )
251  # Old style ImageMapper grid. ZogyTask has its own grid option
252  doSpatiallyVarying = pexConfig.Field(
253  dtype=bool,
254  default=False,
255  doc="Perform A&L decorrelation on a grid across the "
256  "image in order to allow for spatial variations. Zogy does not use this option."
257  )
258  detection = pexConfig.ConfigurableField(
259  target=SourceDetectionTask,
260  doc="Low-threshold detection for final measurement",
261  )
262  measurement = pexConfig.ConfigurableField(
263  target=DipoleFitTask,
264  doc="Enable updated dipole fitting method",
265  )
266  doApCorr = lsst.pex.config.Field(
267  dtype=bool,
268  default=True,
269  doc="Run subtask to apply aperture corrections"
270  )
271  applyApCorr = lsst.pex.config.ConfigurableField(
272  target=ApplyApCorrTask,
273  doc="Subtask to apply aperture corrections"
274  )
275  forcedMeasurement = pexConfig.ConfigurableField(
276  target=ForcedMeasurementTask,
277  doc="Subtask to force photometer PVI at diaSource location.",
278  )
279  getTemplate = pexConfig.ConfigurableField(
280  target=GetCoaddAsTemplateTask,
281  doc="Subtask to retrieve template exposure and sources",
282  )
283  scaleVariance = pexConfig.ConfigurableField(
284  target=ScaleVarianceTask,
285  doc="Subtask to rescale the variance of the template "
286  "to the statistically expected level"
287  )
288  controlStepSize = pexConfig.Field(
289  doc="What step size (every Nth one) to select a control sample from the kernelSources",
290  dtype=int,
291  default=5
292  )
293  controlRandomSeed = pexConfig.Field(
294  doc="Random seed for shuffing the control sample",
295  dtype=int,
296  default=10
297  )
298  register = pexConfig.ConfigurableField(
299  target=RegisterTask,
300  doc="Task to enable image-to-image image registration (warping)",
301  )
302  kernelSourcesFromRef = pexConfig.Field(
303  doc="Select sources to measure kernel from reference catalog if True, template if false",
304  dtype=bool,
305  default=False
306  )
307  templateSipOrder = pexConfig.Field(
308  dtype=int, default=2,
309  doc="Sip Order for fitting the Template Wcs (default is too high, overfitting)"
310  )
311  growFootprint = pexConfig.Field(
312  dtype=int, default=2,
313  doc="Grow positive and negative footprints by this amount before merging"
314  )
315  diaSourceMatchRadius = pexConfig.Field(
316  dtype=float, default=0.5,
317  doc="Match radius (in arcseconds) for DiaSource to Source association"
318  )
319  requiredTemplateFraction = pexConfig.Field(
320  dtype=float, default=0.1,
321  doc="Do not attempt to run task if template covers less than this fraction of pixels."
322  "Setting to 0 will always attempt image subtraction"
323  )
324  doSkySources = pexConfig.Field(
325  dtype=bool,
326  default=False,
327  doc="Generate sky sources?",
328  )
329  skySources = pexConfig.ConfigurableField(
330  target=SkyObjectsTask,
331  doc="Generate sky sources",
332  )
333 
334  def setDefaults(self):
335  # defaults are OK for catalog and diacatalog
336 
337  self.subtract['al'].kernel.name = "AL"
338  self.subtract['al'].kernel.active.fitForBackground = True
339  self.subtract['al'].kernel.active.spatialKernelOrder = 1
340  self.subtract['al'].kernel.active.spatialBgOrder = 2
341 
342  # DiaSource Detection
343  self.detection.thresholdPolarity = "both"
344  self.detection.thresholdValue = 5.0
345  self.detection.reEstimateBackground = False
346  self.detection.thresholdType = "pixel_stdev"
347 
348  # Add filtered flux measurement, the correct measurement for pre-convolved images.
349  # Enable all measurements, regardless of doPreConvolve, as it makes data harvesting easier.
350  # To change that you must modify algorithms.names in the task's applyOverrides method,
351  # after the user has set doPreConvolve.
352  self.measurement.algorithms.names.add('base_PeakLikelihoodFlux')
353  self.measurement.plugins.names |= ['ext_trailedSources_Naive',
354  'base_LocalPhotoCalib',
355  'base_LocalWcs']
356 
357  self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"]
358  self.forcedMeasurement.copyColumns = {
359  "id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"}
360  self.forcedMeasurement.slots.centroid = "base_TransformedCentroid"
361  self.forcedMeasurement.slots.shape = None
362 
363  # For shuffling the control sample
364  random.seed(self.controlRandomSeed)
365 
366  def validate(self):
367  pexConfig.Config.validate(self)
368  if not self.doSubtract and not self.doDetection:
369  raise ValueError("Either doSubtract or doDetection must be enabled.")
370  if self.doMeasurement and not self.doDetection:
371  raise ValueError("Cannot run source measurement without source detection.")
372  if self.doMerge and not self.doDetection:
373  raise ValueError("Cannot run source merging without source detection.")
374  if self.doSkySources and not self.doDetection:
375  raise ValueError("Cannot run sky source creation without source detection.")
376  if self.doUseRegister and not self.doSelectSources:
377  raise ValueError("doUseRegister=True and doSelectSources=False. "
378  "Cannot run RegisterTask without selecting sources.")
379  if hasattr(self.getTemplate, "coaddName"):
380  if self.getTemplate.coaddName != self.coaddName:
381  raise ValueError("Mis-matched coaddName and getTemplate.coaddName in the config.")
382  if self.doScaleDiffimVariance and self.doScaleTemplateVariance:
383  raise ValueError("Scaling the diffim variance and scaling the template variance "
384  "are both set. Please choose one or the other.")
385  # We cannot allow inconsistencies that would lead to None or not available output products
386  if self.subtract.name == 'zogy':
387  if self.doWriteMatchedExp:
388  raise ValueError("doWriteMatchedExp=True Matched exposure is not "
389  "calculated in zogy subtraction.")
390  if self.doAddMetrics:
391  raise ValueError("doAddMetrics=True Kernel metrics does not exist in zogy subtraction.")
392  if self.doDecorrelation:
393  raise ValueError(
394  "doDecorrelation=True The decorrelation afterburner does not exist in zogy subtraction.")
395  if self.doSelectSources:
396  raise ValueError(
397  "doSelectSources=True Selecting sources for PSF matching is not a zogy option.")
398  if self.useGaussianForPreConvolution:
399  raise ValueError(
400  "useGaussianForPreConvolution=True This is an AL subtraction only option.")
401  else:
402  # AL only consistency checks
403  if self.useScoreImageDetection and not self.convolveTemplate:
404  raise ValueError(
405  "convolveTemplate=False and useScoreImageDetection=True "
406  "Pre-convolution and matching of the science image is not a supported operation.")
407  if self.doWriteSubtractedExp and self.useScoreImageDetection:
408  raise ValueError(
409  "doWriteSubtractedExp=True and useScoreImageDetection=True "
410  "Regular difference image is not calculated. "
411  "AL subtraction calculates either the regular difference image or the score image.")
412  if self.doWriteScoreExp and not self.useScoreImageDetection:
413  raise ValueError(
414  "doWriteScoreExp=True and useScoreImageDetection=False "
415  "Score image is not calculated. "
416  "AL subtraction calculates either the regular difference image or the score image.")
417  if self.doAddMetrics and not self.doSubtract:
418  raise ValueError("Subtraction must be enabled for kernel metrics calculation.")
419  if self.useGaussianForPreConvolution and not self.useScoreImageDetection:
420  raise ValueError(
421  "useGaussianForPreConvolution=True and useScoreImageDetection=False "
422  "Gaussian PSF approximation exists only for AL subtraction w/ pre-convolution.")
423 
424 
425 class ImageDifferenceTaskRunner(pipeBase.ButlerInitializedTaskRunner):
426 
427  @staticmethod
428  def getTargetList(parsedCmd, **kwargs):
429  return pipeBase.TaskRunner.getTargetList(parsedCmd, templateIdList=parsedCmd.templateId.idList,
430  **kwargs)
431 
432 
433 class ImageDifferenceTask(pipeBase.CmdLineTask, pipeBase.PipelineTask):
434  """Subtract an image from a template and measure the result
435  """
436  ConfigClass = ImageDifferenceConfig
437  RunnerClass = ImageDifferenceTaskRunner
438  _DefaultName = "imageDifference"
439 
440  def __init__(self, butler=None, **kwargs):
441  """!Construct an ImageDifference Task
442 
443  @param[in] butler Butler object to use in constructing reference object loaders
444  """
445  super().__init__(**kwargs)
446  self.makeSubtask("getTemplate")
447 
448  self.makeSubtask("subtract")
449 
450  if self.config.subtract.name == 'al' and self.config.doDecorrelation:
451  self.makeSubtask("decorrelate")
452 
453  if self.config.doScaleTemplateVariance or self.config.doScaleDiffimVariance:
454  self.makeSubtask("scaleVariance")
455 
456  if self.config.doUseRegister:
457  self.makeSubtask("register")
458  self.schema = afwTable.SourceTable.makeMinimalSchema()
459 
460  if self.config.doSelectSources:
461  self.makeSubtask("sourceSelector")
462  if self.config.kernelSourcesFromRef:
463  self.makeSubtask('refObjLoader', butler=butler)
464  self.makeSubtask("astrometer", refObjLoader=self.refObjLoader)
465 
466  self.algMetadata = dafBase.PropertyList()
467  if self.config.doDetection:
468  self.makeSubtask("detection", schema=self.schema)
469  if self.config.doMeasurement:
470  self.makeSubtask("measurement", schema=self.schema,
471  algMetadata=self.algMetadata)
472  if self.config.doApCorr:
473  self.makeSubtask("applyApCorr", schema=self.measurement.schema)
474  if self.config.doForcedMeasurement:
475  self.schema.addField(
476  "ip_diffim_forced_PsfFlux_instFlux", "D",
477  "Forced PSF flux measured on the direct image.",
478  units="count")
479  self.schema.addField(
480  "ip_diffim_forced_PsfFlux_instFluxErr", "D",
481  "Forced PSF flux error measured on the direct image.",
482  units="count")
483  self.schema.addField(
484  "ip_diffim_forced_PsfFlux_area", "F",
485  "Forced PSF flux effective area of PSF.",
486  units="pixel")
487  self.schema.addField(
488  "ip_diffim_forced_PsfFlux_flag", "Flag",
489  "Forced PSF flux general failure flag.")
490  self.schema.addField(
491  "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag",
492  "Forced PSF flux not enough non-rejected pixels in data to attempt the fit.")
493  self.schema.addField(
494  "ip_diffim_forced_PsfFlux_flag_edge", "Flag",
495  "Forced PSF flux object was too close to the edge of the image to use the full PSF model.")
496  self.makeSubtask("forcedMeasurement", refSchema=self.schema)
497  if self.config.doMatchSources:
498  self.schema.addField("refMatchId", "L", "unique id of reference catalog match")
499  self.schema.addField("srcMatchId", "L", "unique id of source match")
500  if self.config.doSkySources:
501  self.makeSubtask("skySources")
502  self.skySourceKey = self.schema.addField("sky_source", type="Flag", doc="Sky objects.")
503 
504  # initialize InitOutputs
505  self.outputSchema = afwTable.SourceCatalog(self.schema)
506  self.outputSchema.getTable().setMetadata(self.algMetadata)
507 
508  @staticmethod
509  def makeIdFactory(expId, expBits):
510  """Create IdFactory instance for unique 64 bit diaSource id-s.
511 
512  Parameters
513  ----------
514  expId : `int`
515  Exposure id.
516 
517  expBits: `int`
518  Number of used bits in ``expId``.
519 
520  Note
521  ----
522  The diasource id-s consists of the ``expId`` stored fixed in the highest value
523  ``expBits`` of the 64-bit integer plus (bitwise or) a generated sequence number in the
524  low value end of the integer.
525 
526  Returns
527  -------
528  idFactory: `lsst.afw.table.IdFactory`
529  """
530  return ExposureIdInfo(expId, expBits).makeSourceIdFactory()
531 
532  @lsst.utils.inheritDoc(pipeBase.PipelineTask)
533  def runQuantum(self, butlerQC: pipeBase.ButlerQuantumContext,
534  inputRefs: pipeBase.InputQuantizedConnection,
535  outputRefs: pipeBase.OutputQuantizedConnection):
536  inputs = butlerQC.get(inputRefs)
537  self.log.info("Processing %s", butlerQC.quantum.dataId)
538  expId, expBits = butlerQC.quantum.dataId.pack("visit_detector",
539  returnMaxBits=True)
540  idFactory = self.makeIdFactory(expId=expId, expBits=expBits)
541  if self.config.coaddName == 'dcr':
542  templateExposures = inputRefs.dcrCoadds
543  else:
544  templateExposures = inputRefs.coaddExposures
545  templateStruct = self.getTemplate.runQuantum(
546  inputs['exposure'], butlerQC, inputRefs.skyMap, templateExposures
547  )
548 
549  self.checkTemplateIsSufficient(templateStruct.exposure)
550 
551  outputs = self.run(exposure=inputs['exposure'],
552  templateExposure=templateStruct.exposure,
553  idFactory=idFactory)
554  # Consistency with runDataref gen2 handling
555  if outputs.diaSources is None:
556  del outputs.diaSources
557  butlerQC.put(outputs, outputRefs)
558 
559  @timeMethod
560  def runDataRef(self, sensorRef, templateIdList=None):
561  """Subtract an image from a template coadd and measure the result.
562 
563  Data I/O wrapper around `run` using the butler in Gen2.
564 
565  Parameters
566  ----------
567  sensorRef : `lsst.daf.persistence.ButlerDataRef`
568  Sensor-level butler data reference, used for the following data products:
569 
570  Input only:
571  - calexp
572  - psf
573  - ccdExposureId
574  - ccdExposureId_bits
575  - self.config.coaddName + "Coadd_skyMap"
576  - self.config.coaddName + "Coadd"
577  Input or output, depending on config:
578  - self.config.coaddName + "Diff_subtractedExp"
579  Output, depending on config:
580  - self.config.coaddName + "Diff_matchedExp"
581  - self.config.coaddName + "Diff_src"
582 
583  Returns
584  -------
585  results : `lsst.pipe.base.Struct`
586  Returns the Struct by `run`.
587  """
588  subtractedExposureName = self.config.coaddName + "Diff_differenceExp"
589  subtractedExposure = None
590  selectSources = None
591  calexpBackgroundExposure = None
592  self.log.info("Processing %s", sensorRef.dataId)
593 
594  # We make one IdFactory that will be used by both icSrc and src datasets;
595  # I don't know if this is the way we ultimately want to do things, but at least
596  # this ensures the source IDs are fully unique.
597  idFactory = self.makeIdFactory(expId=int(sensorRef.get("ccdExposureId")),
598  expBits=sensorRef.get("ccdExposureId_bits"))
599  if self.config.doAddCalexpBackground:
600  calexpBackgroundExposure = sensorRef.get("calexpBackground")
601 
602  # Retrieve the science image we wish to analyze
603  exposure = sensorRef.get("calexp", immediate=True)
604 
605  # Retrieve the template image
606  template = self.getTemplate.runDataRef(exposure, sensorRef, templateIdList=templateIdList)
607 
608  if sensorRef.datasetExists("src"):
609  self.log.info("Source selection via src product")
610  # Sources already exist; for data release processing
611  selectSources = sensorRef.get("src")
612 
613  if not self.config.doSubtract and self.config.doDetection:
614  # If we don't do subtraction, we need the subtracted exposure from the repo
615  subtractedExposure = sensorRef.get(subtractedExposureName)
616  # Both doSubtract and doDetection cannot be False
617 
618  results = self.run(exposure=exposure,
619  selectSources=selectSources,
620  templateExposure=template.exposure,
621  templateSources=template.sources,
622  idFactory=idFactory,
623  calexpBackgroundExposure=calexpBackgroundExposure,
624  subtractedExposure=subtractedExposure)
625 
626  if self.config.doWriteSources and results.diaSources is not None:
627  sensorRef.put(results.diaSources, self.config.coaddName + "Diff_diaSrc")
628  if self.config.doWriteWarpedExp:
629  sensorRef.put(results.warpedExposure, self.config.coaddName + "Diff_warpedExp")
630  if self.config.doWriteMatchedExp:
631  sensorRef.put(results.matchedExposure, self.config.coaddName + "Diff_matchedExp")
632  if self.config.doAddMetrics and self.config.doSelectSources:
633  sensorRef.put(results.selectSources, self.config.coaddName + "Diff_kernelSrc")
634  if self.config.doWriteSubtractedExp:
635  sensorRef.put(results.subtractedExposure, subtractedExposureName)
636  if self.config.doWriteScoreExp:
637  sensorRef.put(results.scoreExposure, self.config.coaddName + "Diff_scoreExp")
638  return results
639 
640  @timeMethod
641  def run(self, exposure=None, selectSources=None, templateExposure=None, templateSources=None,
642  idFactory=None, calexpBackgroundExposure=None, subtractedExposure=None):
643  """PSF matches, subtract two images and perform detection on the difference image.
644 
645  Parameters
646  ----------
647  exposure : `lsst.afw.image.ExposureF`, optional
648  The science exposure, the minuend in the image subtraction.
649  Can be None only if ``config.doSubtract==False``.
650  selectSources : `lsst.afw.table.SourceCatalog`, optional
651  Identified sources on the science exposure. This catalog is used to
652  select sources in order to perform the AL PSF matching on stamp images
653  around them. The selection steps depend on config options and whether
654  ``templateSources`` and ``matchingSources`` specified.
655  templateExposure : `lsst.afw.image.ExposureF`, optional
656  The template to be subtracted from ``exposure`` in the image subtraction.
657  ``templateExposure`` is modified in place if ``config.doScaleTemplateVariance==True``.
658  The template exposure should cover the same sky area as the science exposure.
659  It is either a stich of patches of a coadd skymap image or a calexp
660  of the same pointing as the science exposure. Can be None only
661  if ``config.doSubtract==False`` and ``subtractedExposure`` is not None.
662  templateSources : `lsst.afw.table.SourceCatalog`, optional
663  Identified sources on the template exposure.
664  idFactory : `lsst.afw.table.IdFactory`
665  Generator object to assign ids to detected sources in the difference image.
666  calexpBackgroundExposure : `lsst.afw.image.ExposureF`, optional
667  Background exposure to be added back to the science exposure
668  if ``config.doAddCalexpBackground==True``
669  subtractedExposure : `lsst.afw.image.ExposureF`, optional
670  If ``config.doSubtract==False`` and ``config.doDetection==True``,
671  performs the post subtraction source detection only on this exposure.
672  Otherwise should be None.
673 
674  Returns
675  -------
676  results : `lsst.pipe.base.Struct`
677  ``subtractedExposure`` : `lsst.afw.image.ExposureF`
678  Difference image.
679  ``scoreExposure`` : `lsst.afw.image.ExposureF` or `None`
680  The zogy score exposure, if calculated.
681  ``matchedExposure`` : `lsst.afw.image.ExposureF`
682  The matched PSF exposure.
683  ``subtractRes`` : `lsst.pipe.base.Struct`
684  The returned result structure of the ImagePsfMatchTask subtask.
685  ``diaSources`` : `lsst.afw.table.SourceCatalog`
686  The catalog of detected sources.
687  ``selectSources`` : `lsst.afw.table.SourceCatalog`
688  The input source catalog with optionally added Qa information.
689 
690  Notes
691  -----
692  The following major steps are included:
693 
694  - warp template coadd to match WCS of image
695  - PSF match image to warped template
696  - subtract image from PSF-matched, warped template
697  - detect sources
698  - measure sources
699 
700  For details about the image subtraction configuration modes
701  see `lsst.ip.diffim`.
702  """
703  subtractRes = None
704  controlSources = None
705  subtractedExposure = None
706  scoreExposure = None
707  diaSources = None
708  kernelSources = None
709  # We'll clone exposure if modified but will still need the original
710  exposureOrig = exposure
711 
712  if self.config.doAddCalexpBackground:
713  mi = exposure.getMaskedImage()
714  mi += calexpBackgroundExposure.getImage()
715 
716  if not exposure.hasPsf():
717  raise pipeBase.TaskError("Exposure has no psf")
718  sciencePsf = exposure.getPsf()
719 
720  if self.config.doSubtract:
721  if self.config.doScaleTemplateVariance:
722  self.log.info("Rescaling template variance")
723  templateVarFactor = self.scaleVariance.run(
724  templateExposure.getMaskedImage())
725  self.log.info("Template variance scaling factor: %.2f", templateVarFactor)
726  self.metadata.add("scaleTemplateVarianceFactor", templateVarFactor)
727  self.metadata.add("psfMatchingAlgorithm", self.config.subtract.name)
728 
729  if self.config.subtract.name == 'zogy':
730  subtractRes = self.subtract.run(exposure, templateExposure, doWarping=True)
731  scoreExposure = subtractRes.scoreExp
732  subtractedExposure = subtractRes.diffExp
733  subtractRes.subtractedExposure = subtractedExposure
734  subtractRes.matchedExposure = None
735 
736  elif self.config.subtract.name == 'al':
737  # compute scienceSigmaOrig: sigma of PSF of science image before pre-convolution
738  # Just need a rough estimate; average positions are fine
739  sciAvgPos = sciencePsf.getAveragePosition()
740  scienceSigmaOrig = sciencePsf.computeShape(sciAvgPos).getDeterminantRadius()
741 
742  templatePsf = templateExposure.getPsf()
743  templateAvgPos = templatePsf.getAveragePosition()
744  templateSigma = templatePsf.computeShape(templateAvgPos).getDeterminantRadius()
745 
746  # if requested, convolve the science exposure with its PSF
747  # (properly, this should be a cross-correlation, but our code does not yet support that)
748  # compute scienceSigmaPost: sigma of science exposure with pre-convolution, if done,
749  # else sigma of original science exposure
750  # TODO: DM-22762 This functional block should be moved into its own method
751  preConvPsf = None
752  if self.config.useScoreImageDetection:
753  self.log.warning("AL likelihood image: pre-convolution of PSF is not implemented.")
754  convControl = afwMath.ConvolutionControl()
755  # cannot convolve in place, so need a new image anyway
756  srcMI = exposure.maskedImage
757  exposure = exposure.clone() # New deep copy
758  srcPsf = sciencePsf
759  if self.config.useGaussianForPreConvolution:
760  self.log.info(
761  "AL likelihood image: Using Gaussian (sigma=%.2f) PSF estimation "
762  "for science image pre-convolution", scienceSigmaOrig)
763  # convolve with a simplified PSF model: a double Gaussian
764  kWidth, kHeight = sciencePsf.getLocalKernel().getDimensions()
765  preConvPsf = SingleGaussianPsf(kWidth, kHeight, scienceSigmaOrig)
766  else:
767  # convolve with science exposure's PSF model
768  self.log.info(
769  "AL likelihood image: Using the science image PSF for pre-convolution.")
770  preConvPsf = srcPsf
771  afwMath.convolve(exposure.maskedImage, srcMI, preConvPsf.getLocalKernel(), convControl)
772  scienceSigmaPost = scienceSigmaOrig*math.sqrt(2)
773  else:
774  scienceSigmaPost = scienceSigmaOrig
775 
776  # If requested, find and select sources from the image
777  # else, AL subtraction will do its own source detection
778  # TODO: DM-22762 This functional block should be moved into its own method
779  if self.config.doSelectSources:
780  if selectSources is None:
781  self.log.warning("Src product does not exist; running detection, measurement,"
782  " selection")
783  # Run own detection and measurement; necessary in nightly processing
784  selectSources = self.subtract.getSelectSources(
785  exposure,
786  sigma=scienceSigmaPost,
787  doSmooth=not self.config.useScoreImageDetection,
788  idFactory=idFactory,
789  )
790 
791  if self.config.doAddMetrics:
792  # Number of basis functions
793 
794  nparam = len(makeKernelBasisList(self.subtract.config.kernel.active,
795  referenceFwhmPix=scienceSigmaPost*FwhmPerSigma,
796  targetFwhmPix=templateSigma*FwhmPerSigma))
797  # Modify the schema of all Sources
798  # DEPRECATED: This is a data dependent (nparam) output product schema
799  # outside the task constructor.
800  # NOTE: The pre-determination of nparam at this point
801  # may be incorrect as the template psf is warped later in
802  # ImagePsfMatchTask.matchExposures()
803  kcQa = KernelCandidateQa(nparam)
804  selectSources = kcQa.addToSchema(selectSources)
805  if self.config.kernelSourcesFromRef:
806  # match exposure sources to reference catalog
807  astromRet = self.astrometer.loadAndMatch(exposure=exposure, sourceCat=selectSources)
808  matches = astromRet.matches
809  elif templateSources:
810  # match exposure sources to template sources
811  mc = afwTable.MatchControl()
812  mc.findOnlyClosest = False
813  matches = afwTable.matchRaDec(templateSources, selectSources, 1.0*geom.arcseconds,
814  mc)
815  else:
816  raise RuntimeError("doSelectSources=True and kernelSourcesFromRef=False,"
817  "but template sources not available. Cannot match science "
818  "sources with template sources. Run process* on data from "
819  "which templates are built.")
820 
821  kernelSources = self.sourceSelector.run(selectSources, exposure=exposure,
822  matches=matches).sourceCat
823  random.shuffle(kernelSources, random.random)
824  controlSources = kernelSources[::self.config.controlStepSize]
825  kernelSources = [k for i, k in enumerate(kernelSources)
826  if i % self.config.controlStepSize]
827 
828  if self.config.doSelectDcrCatalog:
829  redSelector = DiaCatalogSourceSelectorTask(
830  DiaCatalogSourceSelectorConfig(grMin=self.sourceSelector.config.grMax,
831  grMax=99.999))
832  redSources = redSelector.selectStars(exposure, selectSources, matches=matches).starCat
833  controlSources.extend(redSources)
834 
835  blueSelector = DiaCatalogSourceSelectorTask(
836  DiaCatalogSourceSelectorConfig(grMin=-99.999,
837  grMax=self.sourceSelector.config.grMin))
838  blueSources = blueSelector.selectStars(exposure, selectSources,
839  matches=matches).starCat
840  controlSources.extend(blueSources)
841 
842  if self.config.doSelectVariableCatalog:
843  varSelector = DiaCatalogSourceSelectorTask(
844  DiaCatalogSourceSelectorConfig(includeVariable=True))
845  varSources = varSelector.selectStars(exposure, selectSources, matches=matches).starCat
846  controlSources.extend(varSources)
847 
848  self.log.info("Selected %d / %d sources for Psf matching (%d for control sample)",
849  len(kernelSources), len(selectSources), len(controlSources))
850 
851  allresids = {}
852  # TODO: DM-22762 This functional block should be moved into its own method
853  if self.config.doUseRegister:
854  self.log.info("Registering images")
855 
856  if templateSources is None:
857  # Run detection on the template, which is
858  # temporarily background-subtracted
859  # sigma of PSF of template image before warping
860  templateSources = self.subtract.getSelectSources(
861  templateExposure,
862  sigma=templateSigma,
863  doSmooth=True,
864  idFactory=idFactory
865  )
866 
867  # Third step: we need to fit the relative astrometry.
868  #
869  wcsResults = self.fitAstrometry(templateSources, templateExposure, selectSources)
870  warpedExp = self.register.warpExposure(templateExposure, wcsResults.wcs,
871  exposure.getWcs(), exposure.getBBox())
872  templateExposure = warpedExp
873 
874  # Create debugging outputs on the astrometric
875  # residuals as a function of position. Persistence
876  # not yet implemented; expected on (I believe) #2636.
877  if self.config.doDebugRegister:
878  # Grab matches to reference catalog
879  srcToMatch = {x.second.getId(): x.first for x in matches}
880 
881  refCoordKey = wcsResults.matches[0].first.getTable().getCoordKey()
882  inCentroidKey = wcsResults.matches[0].second.getTable().getCentroidSlot().getMeasKey()
883  sids = [m.first.getId() for m in wcsResults.matches]
884  positions = [m.first.get(refCoordKey) for m in wcsResults.matches]
885  residuals = [m.first.get(refCoordKey).getOffsetFrom(wcsResults.wcs.pixelToSky(
886  m.second.get(inCentroidKey))) for m in wcsResults.matches]
887  allresids = dict(zip(sids, zip(positions, residuals)))
888 
889  cresiduals = [m.first.get(refCoordKey).getTangentPlaneOffset(
890  wcsResults.wcs.pixelToSky(
891  m.second.get(inCentroidKey))) for m in wcsResults.matches]
892  colors = numpy.array([-2.5*numpy.log10(srcToMatch[x].get("g"))
893  + 2.5*numpy.log10(srcToMatch[x].get("r"))
894  for x in sids if x in srcToMatch.keys()])
895  dlong = numpy.array([r[0].asArcseconds() for s, r in zip(sids, cresiduals)
896  if s in srcToMatch.keys()])
897  dlat = numpy.array([r[1].asArcseconds() for s, r in zip(sids, cresiduals)
898  if s in srcToMatch.keys()])
899  idx1 = numpy.where(colors < self.sourceSelector.config.grMin)
900  idx2 = numpy.where((colors >= self.sourceSelector.config.grMin)
901  & (colors <= self.sourceSelector.config.grMax))
902  idx3 = numpy.where(colors > self.sourceSelector.config.grMax)
903  rms1Long = IqrToSigma*(
904  (numpy.percentile(dlong[idx1], 75) - numpy.percentile(dlong[idx1], 25)))
905  rms1Lat = IqrToSigma*(numpy.percentile(dlat[idx1], 75)
906  - numpy.percentile(dlat[idx1], 25))
907  rms2Long = IqrToSigma*(
908  (numpy.percentile(dlong[idx2], 75) - numpy.percentile(dlong[idx2], 25)))
909  rms2Lat = IqrToSigma*(numpy.percentile(dlat[idx2], 75)
910  - numpy.percentile(dlat[idx2], 25))
911  rms3Long = IqrToSigma*(
912  (numpy.percentile(dlong[idx3], 75) - numpy.percentile(dlong[idx3], 25)))
913  rms3Lat = IqrToSigma*(numpy.percentile(dlat[idx3], 75)
914  - numpy.percentile(dlat[idx3], 25))
915  self.log.info("Blue star offsets'': %.3f %.3f, %.3f %.3f",
916  numpy.median(dlong[idx1]), rms1Long,
917  numpy.median(dlat[idx1]), rms1Lat)
918  self.log.info("Green star offsets'': %.3f %.3f, %.3f %.3f",
919  numpy.median(dlong[idx2]), rms2Long,
920  numpy.median(dlat[idx2]), rms2Lat)
921  self.log.info("Red star offsets'': %.3f %.3f, %.3f %.3f",
922  numpy.median(dlong[idx3]), rms3Long,
923  numpy.median(dlat[idx3]), rms3Lat)
924 
925  self.metadata.add("RegisterBlueLongOffsetMedian", numpy.median(dlong[idx1]))
926  self.metadata.add("RegisterGreenLongOffsetMedian", numpy.median(dlong[idx2]))
927  self.metadata.add("RegisterRedLongOffsetMedian", numpy.median(dlong[idx3]))
928  self.metadata.add("RegisterBlueLongOffsetStd", rms1Long)
929  self.metadata.add("RegisterGreenLongOffsetStd", rms2Long)
930  self.metadata.add("RegisterRedLongOffsetStd", rms3Long)
931 
932  self.metadata.add("RegisterBlueLatOffsetMedian", numpy.median(dlat[idx1]))
933  self.metadata.add("RegisterGreenLatOffsetMedian", numpy.median(dlat[idx2]))
934  self.metadata.add("RegisterRedLatOffsetMedian", numpy.median(dlat[idx3]))
935  self.metadata.add("RegisterBlueLatOffsetStd", rms1Lat)
936  self.metadata.add("RegisterGreenLatOffsetStd", rms2Lat)
937  self.metadata.add("RegisterRedLatOffsetStd", rms3Lat)
938 
939  # warp template exposure to match exposure,
940  # PSF match template exposure to exposure,
941  # then return the difference
942 
943  # Return warped template... Construct sourceKernelCand list after subtract
944  self.log.info("Subtracting images")
945  subtractRes = self.subtract.subtractExposures(
946  templateExposure=templateExposure,
947  scienceExposure=exposure,
948  candidateList=kernelSources,
949  convolveTemplate=self.config.convolveTemplate,
950  doWarping=not self.config.doUseRegister
951  )
952  if self.config.useScoreImageDetection:
953  scoreExposure = subtractRes.subtractedExposure
954  else:
955  subtractedExposure = subtractRes.subtractedExposure
956 
957  if self.config.doDetection:
958  self.log.info("Computing diffim PSF")
959 
960  # Get Psf from the appropriate input image if it doesn't exist
961  if subtractedExposure is not None and not subtractedExposure.hasPsf():
962  if self.config.convolveTemplate:
963  subtractedExposure.setPsf(exposure.getPsf())
964  else:
965  subtractedExposure.setPsf(templateExposure.getPsf())
966 
967  # If doSubtract is False, then subtractedExposure was fetched from disk (above),
968  # thus it may have already been decorrelated. Thus, we do not decorrelate if
969  # doSubtract is False.
970 
971  # NOTE: At this point doSubtract == True
972  if self.config.doDecorrelation and self.config.doSubtract:
973  preConvKernel = None
974  if self.config.useGaussianForPreConvolution:
975  preConvKernel = preConvPsf.getLocalKernel()
976  if self.config.useScoreImageDetection:
977  scoreExposure = self.decorrelate.run(exposureOrig, subtractRes.warpedExposure,
978  scoreExposure,
979  subtractRes.psfMatchingKernel,
980  spatiallyVarying=self.config.doSpatiallyVarying,
981  preConvKernel=preConvKernel,
982  templateMatched=True,
983  preConvMode=True).correctedExposure
984  # Note that the subtracted exposure is always decorrelated,
985  # even if the score image is used for detection
986  subtractedExposure = self.decorrelate.run(exposureOrig, subtractRes.warpedExposure,
987  subtractedExposure,
988  subtractRes.psfMatchingKernel,
989  spatiallyVarying=self.config.doSpatiallyVarying,
990  preConvKernel=None,
991  templateMatched=self.config.convolveTemplate,
992  preConvMode=False).correctedExposure
993  # END (if subtractAlgorithm == 'AL')
994  # END (if self.config.doSubtract)
995  if self.config.doDetection:
996  self.log.info("Running diaSource detection")
997 
998  # subtractedExposure - reserved for task return value
999  # in zogy, it is always the proper difference image
1000  # in AL, it may be (yet) pre-convolved and/or decorrelated
1001  #
1002  # detectionExposure - controls which exposure to use for detection
1003  # in-place modifications will appear in task return
1004  if self.config.useScoreImageDetection:
1005  # zogy with score image detection enabled
1006  self.log.info("Detection, diffim rescaling and measurements are "
1007  "on AL likelihood or Zogy score image.")
1008  detectionExposure = scoreExposure
1009  else:
1010  # AL or zogy with no score image detection
1011  detectionExposure = subtractedExposure
1012 
1013  # Rescale difference image variance plane
1014  if self.config.doScaleDiffimVariance:
1015  self.log.info("Rescaling diffim variance")
1016  diffimVarFactor = self.scaleVariance.run(detectionExposure.getMaskedImage())
1017  self.log.info("Diffim variance scaling factor: %.2f", diffimVarFactor)
1018  self.metadata.add("scaleDiffimVarianceFactor", diffimVarFactor)
1019 
1020  # Erase existing detection mask planes
1021  mask = detectionExposure.getMaskedImage().getMask()
1022  mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
1023 
1024  table = afwTable.SourceTable.make(self.schema, idFactory)
1025  table.setMetadata(self.algMetadata)
1026  results = self.detection.run(
1027  table=table,
1028  exposure=detectionExposure,
1029  doSmooth=not self.config.useScoreImageDetection
1030  )
1031 
1032  if self.config.doMerge:
1033  fpSet = results.fpSets.positive
1034  fpSet.merge(results.fpSets.negative, self.config.growFootprint,
1035  self.config.growFootprint, False)
1036  diaSources = afwTable.SourceCatalog(table)
1037  fpSet.makeSources(diaSources)
1038  self.log.info("Merging detections into %d sources", len(diaSources))
1039  else:
1040  diaSources = results.sources
1041  # Inject skySources before measurement.
1042  if self.config.doSkySources:
1043  skySourceFootprints = self.skySources.run(
1044  mask=detectionExposure.mask,
1045  seed=detectionExposure.info.id)
1046  if skySourceFootprints:
1047  for foot in skySourceFootprints:
1048  s = diaSources.addNew()
1049  s.setFootprint(foot)
1050  s.set(self.skySourceKey, True)
1051 
1052  if self.config.doMeasurement:
1053  newDipoleFitting = self.config.doDipoleFitting
1054  self.log.info("Running diaSource measurement: newDipoleFitting=%r", newDipoleFitting)
1055  if not newDipoleFitting:
1056  # Just fit dipole in diffim
1057  self.measurement.run(diaSources, detectionExposure)
1058  else:
1059  # Use (matched) template and science image (if avail.) to constrain dipole fitting
1060  if self.config.doSubtract and 'matchedExposure' in subtractRes.getDict():
1061  self.measurement.run(diaSources, detectionExposure, exposure,
1062  subtractRes.matchedExposure)
1063  else:
1064  self.measurement.run(diaSources, detectionExposure, exposure)
1065  if self.config.doApCorr:
1066  self.applyApCorr.run(
1067  catalog=diaSources,
1068  apCorrMap=detectionExposure.getInfo().getApCorrMap()
1069  )
1070 
1071  if self.config.doForcedMeasurement:
1072  # Run forced psf photometry on the PVI at the diaSource locations.
1073  # Copy the measured flux and error into the diaSource.
1074  forcedSources = self.forcedMeasurement.generateMeasCat(
1075  exposure, diaSources, detectionExposure.getWcs())
1076  self.forcedMeasurement.run(forcedSources, exposure, diaSources, detectionExposure.getWcs())
1077  mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
1078  mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFlux")[0],
1079  "ip_diffim_forced_PsfFlux_instFlux", True)
1080  mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFluxErr")[0],
1081  "ip_diffim_forced_PsfFlux_instFluxErr", True)
1082  mapper.addMapping(forcedSources.schema.find("base_PsfFlux_area")[0],
1083  "ip_diffim_forced_PsfFlux_area", True)
1084  mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag")[0],
1085  "ip_diffim_forced_PsfFlux_flag", True)
1086  mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_noGoodPixels")[0],
1087  "ip_diffim_forced_PsfFlux_flag_noGoodPixels", True)
1088  mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_edge")[0],
1089  "ip_diffim_forced_PsfFlux_flag_edge", True)
1090  for diaSource, forcedSource in zip(diaSources, forcedSources):
1091  diaSource.assign(forcedSource, mapper)
1092 
1093  # Match with the calexp sources if possible
1094  if self.config.doMatchSources:
1095  if selectSources is not None:
1096  # Create key,val pair where key=diaSourceId and val=sourceId
1097  matchRadAsec = self.config.diaSourceMatchRadius
1098  matchRadPixel = matchRadAsec/exposure.getWcs().getPixelScale().asArcseconds()
1099 
1100  srcMatches = afwTable.matchXy(selectSources, diaSources, matchRadPixel)
1101  srcMatchDict = dict([(srcMatch.second.getId(), srcMatch.first.getId()) for
1102  srcMatch in srcMatches])
1103  self.log.info("Matched %d / %d diaSources to sources",
1104  len(srcMatchDict), len(diaSources))
1105  else:
1106  self.log.warning("Src product does not exist; cannot match with diaSources")
1107  srcMatchDict = {}
1108 
1109  # Create key,val pair where key=diaSourceId and val=refId
1110  refAstromConfig = AstrometryConfig()
1111  refAstromConfig.matcher.maxMatchDistArcSec = matchRadAsec
1112  refAstrometer = AstrometryTask(refAstromConfig)
1113  astromRet = refAstrometer.run(exposure=exposure, sourceCat=diaSources)
1114  refMatches = astromRet.matches
1115  if refMatches is None:
1116  self.log.warning("No diaSource matches with reference catalog")
1117  refMatchDict = {}
1118  else:
1119  self.log.info("Matched %d / %d diaSources to reference catalog",
1120  len(refMatches), len(diaSources))
1121  refMatchDict = dict([(refMatch.second.getId(), refMatch.first.getId()) for
1122  refMatch in refMatches])
1123 
1124  # Assign source Ids
1125  for diaSource in diaSources:
1126  sid = diaSource.getId()
1127  if sid in srcMatchDict:
1128  diaSource.set("srcMatchId", srcMatchDict[sid])
1129  if sid in refMatchDict:
1130  diaSource.set("refMatchId", refMatchDict[sid])
1131 
1132  if self.config.doAddMetrics and self.config.doSelectSources:
1133  self.log.info("Evaluating metrics and control sample")
1134 
1135  kernelCandList = []
1136  for cell in subtractRes.kernelCellSet.getCellList():
1137  for cand in cell.begin(False): # include bad candidates
1138  kernelCandList.append(cand)
1139 
1140  # Get basis list to build control sample kernels
1141  basisList = kernelCandList[0].getKernel(KernelCandidateF.ORIG).getKernelList()
1142  nparam = len(kernelCandList[0].getKernel(KernelCandidateF.ORIG).getKernelParameters())
1143 
1144  controlCandList = (
1145  diffimTools.sourceTableToCandidateList(controlSources,
1146  subtractRes.warpedExposure, exposure,
1147  self.config.subtract.kernel.active,
1148  self.config.subtract.kernel.active.detectionConfig,
1149  self.log, doBuild=True, basisList=basisList))
1150 
1151  KernelCandidateQa.apply(kernelCandList, subtractRes.psfMatchingKernel,
1152  subtractRes.backgroundModel, dof=nparam)
1153  KernelCandidateQa.apply(controlCandList, subtractRes.psfMatchingKernel,
1154  subtractRes.backgroundModel)
1155 
1156  if self.config.doDetection:
1157  KernelCandidateQa.aggregate(selectSources, self.metadata, allresids, diaSources)
1158  else:
1159  KernelCandidateQa.aggregate(selectSources, self.metadata, allresids)
1160 
1161  self.runDebug(exposure, subtractRes, selectSources, kernelSources, diaSources)
1162  return pipeBase.Struct(
1163  subtractedExposure=subtractedExposure,
1164  scoreExposure=scoreExposure,
1165  warpedExposure=subtractRes.warpedExposure,
1166  matchedExposure=subtractRes.matchedExposure,
1167  subtractRes=subtractRes,
1168  diaSources=diaSources,
1169  selectSources=selectSources
1170  )
1171 
1172  def fitAstrometry(self, templateSources, templateExposure, selectSources):
1173  """Fit the relative astrometry between templateSources and selectSources
1174 
1175  Todo
1176  ----
1177 
1178  Remove this method. It originally fit a new WCS to the template before calling register.run
1179  because our TAN-SIP fitter behaved badly for points far from CRPIX, but that's been fixed.
1180  It remains because a subtask overrides it.
1181  """
1182  results = self.register.run(templateSources, templateExposure.getWcs(),
1183  templateExposure.getBBox(), selectSources)
1184  return results
1185 
1186  def runDebug(self, exposure, subtractRes, selectSources, kernelSources, diaSources):
1187  """Make debug plots and displays.
1188 
1189  Todo
1190  ----
1191  Test and update for current debug display and slot names
1192  """
1193  import lsstDebug
1194  display = lsstDebug.Info(__name__).display
1195  showSubtracted = lsstDebug.Info(__name__).showSubtracted
1196  showPixelResiduals = lsstDebug.Info(__name__).showPixelResiduals
1197  showDiaSources = lsstDebug.Info(__name__).showDiaSources
1198  showDipoles = lsstDebug.Info(__name__).showDipoles
1199  maskTransparency = lsstDebug.Info(__name__).maskTransparency
1200  if display:
1201  disp = afwDisplay.getDisplay(frame=lsstDebug.frame)
1202  if not maskTransparency:
1203  maskTransparency = 0
1204  disp.setMaskTransparency(maskTransparency)
1205 
1206  if display and showSubtracted:
1207  disp.mtv(subtractRes.subtractedExposure, title="Subtracted image")
1208  mi = subtractRes.subtractedExposure.getMaskedImage()
1209  x0, y0 = mi.getX0(), mi.getY0()
1210  with disp.Buffering():
1211  for s in diaSources:
1212  x, y = s.getX() - x0, s.getY() - y0
1213  ctype = "red" if s.get("flags_negative") else "yellow"
1214  if (s.get("base_PixelFlags_flag_interpolatedCenter")
1215  or s.get("base_PixelFlags_flag_saturatedCenter")
1216  or s.get("base_PixelFlags_flag_crCenter")):
1217  ptype = "x"
1218  elif (s.get("base_PixelFlags_flag_interpolated")
1219  or s.get("base_PixelFlags_flag_saturated")
1220  or s.get("base_PixelFlags_flag_cr")):
1221  ptype = "+"
1222  else:
1223  ptype = "o"
1224  disp.dot(ptype, x, y, size=4, ctype=ctype)
1225  lsstDebug.frame += 1
1226 
1227  if display and showPixelResiduals and selectSources:
1228  nonKernelSources = []
1229  for source in selectSources:
1230  if source not in kernelSources:
1231  nonKernelSources.append(source)
1232 
1233  diUtils.plotPixelResiduals(exposure,
1234  subtractRes.warpedExposure,
1235  subtractRes.subtractedExposure,
1236  subtractRes.kernelCellSet,
1237  subtractRes.psfMatchingKernel,
1238  subtractRes.backgroundModel,
1239  nonKernelSources,
1240  self.subtract.config.kernel.active.detectionConfig,
1241  origVariance=False)
1242  diUtils.plotPixelResiduals(exposure,
1243  subtractRes.warpedExposure,
1244  subtractRes.subtractedExposure,
1245  subtractRes.kernelCellSet,
1246  subtractRes.psfMatchingKernel,
1247  subtractRes.backgroundModel,
1248  nonKernelSources,
1249  self.subtract.config.kernel.active.detectionConfig,
1250  origVariance=True)
1251  if display and showDiaSources:
1252  flagChecker = SourceFlagChecker(diaSources)
1253  isFlagged = [flagChecker(x) for x in diaSources]
1254  isDipole = [x.get("ip_diffim_ClassificationDipole_value") for x in diaSources]
1255  diUtils.showDiaSources(diaSources, subtractRes.subtractedExposure, isFlagged, isDipole,
1256  frame=lsstDebug.frame)
1257  lsstDebug.frame += 1
1258 
1259  if display and showDipoles:
1260  DipoleAnalysis().displayDipoles(subtractRes.subtractedExposure, diaSources,
1261  frame=lsstDebug.frame)
1262  lsstDebug.frame += 1
1263 
1264  def checkTemplateIsSufficient(self, templateExposure):
1265  """Raise NoWorkFound if template coverage < requiredTemplateFraction
1266 
1267  Parameters
1268  ----------
1269  templateExposure : `lsst.afw.image.ExposureF`
1270  The template exposure to check
1271 
1272  Raises
1273  ------
1274  NoWorkFound
1275  Raised if fraction of good pixels, defined as not having NO_DATA
1276  set, is less then the configured requiredTemplateFraction
1277  """
1278  # Count the number of pixels with the NO_DATA mask bit set
1279  # counting NaN pixels is insufficient because pixels without data are often intepolated over)
1280  pixNoData = numpy.count_nonzero(templateExposure.mask.array
1281  & templateExposure.mask.getPlaneBitMask('NO_DATA'))
1282  pixGood = templateExposure.getBBox().getArea() - pixNoData
1283  self.log.info("template has %d good pixels (%.1f%%)", pixGood,
1284  100*pixGood/templateExposure.getBBox().getArea())
1285 
1286  if pixGood/templateExposure.getBBox().getArea() < self.config.requiredTemplateFraction:
1287  message = ("Insufficient Template Coverage. (%.1f%% < %.1f%%) Not attempting subtraction. "
1288  "To force subtraction, set config requiredTemplateFraction=0." % (
1289  100*pixGood/templateExposure.getBBox().getArea(),
1290  100*self.config.requiredTemplateFraction))
1291  raise pipeBase.NoWorkFound(message)
1292 
1293  def _getConfigName(self):
1294  """Return the name of the config dataset
1295  """
1296  return "%sDiff_config" % (self.config.coaddName,)
1297 
1298  def _getMetadataName(self):
1299  """Return the name of the metadata dataset
1300  """
1301  return "%sDiff_metadata" % (self.config.coaddName,)
1302 
1303  def getSchemaCatalogs(self):
1304  """Return a dict of empty catalogs for each catalog dataset produced by this task."""
1305  return {self.config.coaddName + "Diff_diaSrc": self.outputSchema}
1306 
1307  @classmethod
1308  def _makeArgumentParser(cls):
1309  """Create an argument parser
1310  """
1311  parser = pipeBase.ArgumentParser(name=cls._DefaultName)
1312  parser.add_id_argument("--id", "calexp", help="data ID, e.g. --id visit=12345 ccd=1,2")
1313  parser.add_id_argument("--templateId", "calexp", doMakeDataRefList=True,
1314  help="Template data ID in case of calexp template,"
1315  " e.g. --templateId visit=6789")
1316  return parser
1317 
1318 
1319 class ImageDifferenceFromTemplateConnections(ImageDifferenceTaskConnections,
1320  defaultTemplates={"coaddName": "goodSeeing"}
1321  ):
1322  inputTemplate = pipeBase.connectionTypes.Input(
1323  doc=("Warped template produced by GetMultiTractCoaddTemplate"),
1324  dimensions=("instrument", "visit", "detector"),
1325  storageClass="ExposureF",
1326  name="{fakesType}{coaddName}Diff_templateExp{warpTypeSuffix}",
1327  )
1328 
1329  def __init__(self, *, config=None):
1330  super().__init__(config=config)
1331  # ImageDifferenceConnections will have removed one of these.
1332  # Make sure they're both gone, because no coadds are needed.
1333  if "coaddExposures" in self.inputs:
1334  self.inputs.remove("coaddExposures")
1335  if "dcrCoadds" in self.inputs:
1336  self.inputs.remove("dcrCoadds")
1337 
1338 
1339 class ImageDifferenceFromTemplateConfig(ImageDifferenceConfig,
1340  pipelineConnections=ImageDifferenceFromTemplateConnections):
1341  pass
1342 
1343 
1344 class ImageDifferenceFromTemplateTask(ImageDifferenceTask):
1345  ConfigClass = ImageDifferenceFromTemplateConfig
1346  _DefaultName = "imageDifference"
1347 
1348  @lsst.utils.inheritDoc(pipeBase.PipelineTask)
1349  def runQuantum(self, butlerQC, inputRefs, outputRefs):
1350  inputs = butlerQC.get(inputRefs)
1351  self.log.info("Processing %s", butlerQC.quantum.dataId)
1352  self.checkTemplateIsSufficient(inputs['inputTemplate'])
1353  expId, expBits = butlerQC.quantum.dataId.pack("visit_detector",
1354  returnMaxBits=True)
1355  idFactory = self.makeIdFactory(expId=expId, expBits=expBits)
1356 
1357  outputs = self.run(exposure=inputs['exposure'],
1358  templateExposure=inputs['inputTemplate'],
1359  idFactory=idFactory)
1360 
1361  # Consistency with runDataref gen2 handling
1362  if outputs.diaSources is None:
1363  del outputs.diaSources
1364  butlerQC.put(outputs, outputRefs)
1365 
1366 
1367 class Winter2013ImageDifferenceConfig(ImageDifferenceConfig):
1368  winter2013WcsShift = pexConfig.Field(dtype=float, default=0.0,
1369  doc="Shift stars going into RegisterTask by this amount")
1370  winter2013WcsRms = pexConfig.Field(dtype=float, default=0.0,
1371  doc="Perturb stars going into RegisterTask by this amount")
1372 
1373  def setDefaults(self):
1374  ImageDifferenceConfig.setDefaults(self)
1375  self.getTemplate.retarget(GetCalexpAsTemplateTask)
1376 
1377 
1378 class Winter2013ImageDifferenceTask(ImageDifferenceTask):
1379  """!Image difference Task used in the Winter 2013 data challege.
1380  Enables testing the effects of registration shifts and scatter.
1381 
1382  For use with winter 2013 simulated images:
1383  Use --templateId visit=88868666 for sparse data
1384  --templateId visit=22222200 for dense data (g)
1385  --templateId visit=11111100 for dense data (i)
1386  """
1387  ConfigClass = Winter2013ImageDifferenceConfig
1388  _DefaultName = "winter2013ImageDifference"
1389 
1390  def __init__(self, **kwargs):
1391  ImageDifferenceTask.__init__(self, **kwargs)
1392 
1393  def fitAstrometry(self, templateSources, templateExposure, selectSources):
1394  """Fit the relative astrometry between templateSources and selectSources"""
1395  if self.config.winter2013WcsShift > 0.0:
1396  offset = geom.Extent2D(self.config.winter2013WcsShift,
1397  self.config.winter2013WcsShift)
1398  cKey = templateSources[0].getTable().getCentroidSlot().getMeasKey()
1399  for source in templateSources:
1400  centroid = source.get(cKey)
1401  source.set(cKey, centroid + offset)
1402  elif self.config.winter2013WcsRms > 0.0:
1403  cKey = templateSources[0].getTable().getCentroidSlot().getMeasKey()
1404  for source in templateSources:
1405  offset = geom.Extent2D(self.config.winter2013WcsRms*numpy.random.normal(),
1406  self.config.winter2013WcsRms*numpy.random.normal())
1407  centroid = source.get(cKey)
1408  source.set(cKey, centroid + offset)
1409 
1410  results = self.register.run(templateSources, templateExposure.getWcs(),
1411  templateExposure.getBBox(), selectSources)
1412  return results
Parameters to control convolution.
Definition: ConvolveImage.h:50
Pass parameters to algorithms that match list of sources.
Definition: Match.h:45
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
Represent a PSF as a circularly symmetrical Gaussian.
def runQuantum(self, butlerQC, inputRefs, outputRefs)
Image difference Task used in the Winter 2013 data challege.
def fitAstrometry(self, templateSources, templateExposure, selectSources)
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
int warpExposure(DestExposureT &destExposure, SrcExposureT const &srcExposure, WarpingControl const &control, typename DestExposureT::MaskedImageT::SinglePixel padValue=lsst::afw::math::edgePixel< typename DestExposureT::MaskedImageT >(typename lsst::afw::image::detail::image_traits< typename DestExposureT::MaskedImageT >::image_category()))
Warp (remap) one exposure to another.
SourceMatchVector matchXy(SourceCatalog const &cat1, SourceCatalog const &cat2, double radius, MatchControl const &mc=MatchControl())
Compute all tuples (s1,s2,d) where s1 belings to cat1, s2 belongs to cat2 and d, the distance between...
Definition: Match.cc:305
std::vector< Match< typename Cat1::Record, typename Cat2::Record > > matchRaDec(Cat1 const &cat1, Cat2 const &cat2, lsst::geom::Angle radius, MatchControl const &mc=MatchControl())
Compute all tuples (s1,s2,d) where s1 belings to cat1, s2 belongs to cat2 and d, the distance between...
Definition: Match.cc:158
def run(self, coaddExposures, bbox, wcs)
Definition: getTemplate.py:603
def makeKernelBasisList(config, targetFwhmPix=None, referenceFwhmPix=None, basisDegGauss=None, basisSigmaGauss=None, metadata=None)
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.