LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
getTemplate.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2016 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 
23 import numpy as np
24 
25 import lsst.afw.image as afwImage
26 import lsst.geom as geom
27 import lsst.afw.geom as afwGeom
28 import lsst.afw.table as afwTable
29 import lsst.afw.math as afwMath
30 import lsst.pex.config as pexConfig
31 import lsst.pipe.base as pipeBase
32 from lsst.skymap import BaseSkyMap
33 from lsst.daf.butler import DeferredDatasetHandle
34 from lsst.ip.diffim.dcrModel import DcrModel
35 from lsst.meas.algorithms import CoaddPsf, WarpedPsf, CoaddPsfConfig
36 
37 __all__ = ["GetCoaddAsTemplateTask", "GetCoaddAsTemplateConfig",
38  "GetCalexpAsTemplateTask", "GetCalexpAsTemplateConfig",
39  "GetMultiTractCoaddTemplateTask", "GetMultiTractCoaddTemplateConfig"]
40 
41 
42 class GetCoaddAsTemplateConfig(pexConfig.Config):
43  templateBorderSize = pexConfig.Field(
44  dtype=int,
45  default=10,
46  doc="Number of pixels to grow the requested template image to account for warping"
47  )
48  coaddName = pexConfig.Field(
49  doc="coadd name: typically one of 'deep', 'goodSeeing', or 'dcr'",
50  dtype=str,
51  default="deep",
52  )
53  numSubfilters = pexConfig.Field(
54  doc="Number of subfilters in the DcrCoadd. Used only if ``coaddName``='dcr'",
55  dtype=int,
56  default=3,
57  )
58  effectiveWavelength = pexConfig.Field(
59  doc="Effective wavelength of the filter. Used only if ``coaddName``='dcr'",
60  optional=True,
61  dtype=float,
62  )
63  bandwidth = pexConfig.Field(
64  doc="Bandwidth of the physical filter. Used only if ``coaddName``='dcr'",
65  optional=True,
66  dtype=float,
67  )
68  warpType = pexConfig.Field(
69  doc="Warp type of the coadd template: one of 'direct' or 'psfMatched'",
70  dtype=str,
71  default="direct",
72  )
73 
74  def validate(self):
75  if self.coaddNamecoaddName == 'dcr':
76  if self.effectiveWavelengtheffectiveWavelength is None or self.bandwidthbandwidth is None:
77  raise ValueError("The effective wavelength and bandwidth of the physical filter "
78  "must be set in the getTemplate config for DCR coadds. "
79  "Required until transmission curves are used in DM-13668.")
80 
81 
82 class GetCoaddAsTemplateTask(pipeBase.Task):
83  """Subtask to retrieve coadd for use as an image difference template.
84 
85  This is the default getTemplate Task to be run as a subtask by
86  ``pipe.tasks.ImageDifferenceTask``. The main methods are ``run()`` and
87  ``runGen3()``.
88 
89  Notes
90  -----
91  From the given skymap, the closest tract is selected; multiple tracts are
92  not supported. The assembled template inherits the WCS of the selected
93  skymap tract and the resolution of the template exposures. Overlapping box
94  regions of the input template patches are pixel by pixel copied into the
95  assembled template image. There is no warping or pixel resampling.
96 
97  Pixels with no overlap of any available input patches are set to ``nan`` value
98  and ``NO_DATA`` flagged.
99  """
100 
101  ConfigClass = GetCoaddAsTemplateConfig
102  _DefaultName = "GetCoaddAsTemplateTask"
103 
104  def runDataRef(self, exposure, sensorRef, templateIdList=None):
105  """Gen2 task entry point. Retrieve and mosaic a template coadd exposure
106  that overlaps the science exposure.
107 
108  Parameters
109  ----------
110  exposure: `lsst.afw.image.Exposure`
111  an exposure for which to generate an overlapping template
112  sensorRef : TYPE
113  a Butler data reference that can be used to obtain coadd data
114  templateIdList : TYPE, optional
115  list of data ids, unused here, in the case of coadd template
116 
117  Returns
118  -------
119  result : `lsst.pipe.base.Struct`
120  - ``exposure`` : `lsst.afw.image.ExposureF`
121  a template coadd exposure assembled out of patches
122  - ``sources`` : None for this subtask
123  """
124  skyMap = sensorRef.get(datasetType=self.config.coaddName + "Coadd_skyMap")
125  tractInfo, patchList, skyCorners = self.getOverlapPatchListgetOverlapPatchList(exposure, skyMap)
126 
127  availableCoaddRefs = dict()
128  for patchInfo in patchList:
129  patchNumber = tractInfo.getSequentialPatchIndex(patchInfo)
130  patchArgDict = dict(
131  datasetType=self.getCoaddDatasetNamegetCoaddDatasetName() + "_sub",
132  bbox=patchInfo.getOuterBBox(),
133  tract=tractInfo.getId(),
134  patch="%s,%s" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),
135  subfilter=0,
136  numSubfilters=self.config.numSubfilters,
137  )
138 
139  if sensorRef.datasetExists(**patchArgDict):
140  self.log.info("Reading patch %s", patchArgDict)
141  availableCoaddRefs[patchNumber] = patchArgDict
142 
143  templateExposure = self.runrun(
144  tractInfo, patchList, skyCorners, availableCoaddRefs,
145  sensorRef=sensorRef, visitInfo=exposure.getInfo().getVisitInfo()
146  )
147  return pipeBase.Struct(exposure=templateExposure, sources=None)
148 
149  def runQuantum(self, exposure, butlerQC, skyMapRef, coaddExposureRefs):
150  """Gen3 task entry point. Retrieve and mosaic a template coadd exposure
151  that overlaps the science exposure.
152 
153  Parameters
154  ----------
155  exposure : `lsst.afw.image.Exposure`
156  The science exposure to define the sky region of the template coadd.
157  butlerQC : `lsst.pipe.base.ButlerQuantumContext`
158  Butler like object that supports getting data by DatasetRef.
159  skyMapRef : `lsst.daf.butler.DatasetRef`
160  Reference to SkyMap object that corresponds to the template coadd.
161  coaddExposureRefs : iterable of `lsst.daf.butler.DeferredDatasetRef`
162  Iterable of references to the available template coadd patches.
163 
164  Returns
165  -------
166  result : `lsst.pipe.base.Struct`
167  - ``exposure`` : `lsst.afw.image.ExposureF`
168  a template coadd exposure assembled out of patches
169  - ``sources`` : `None` for this subtask
170  """
171  skyMap = butlerQC.get(skyMapRef)
172  coaddExposureRefs = butlerQC.get(coaddExposureRefs)
173  tracts = [ref.dataId['tract'] for ref in coaddExposureRefs]
174  if tracts.count(tracts[0]) == len(tracts):
175  tractInfo = skyMap[tracts[0]]
176  else:
177  raise RuntimeError("Templates constructed from multiple Tracts not supported by this task. "
178  "Use GetMultiTractCoaddTemplateTask instead.")
179 
180  detectorBBox = exposure.getBBox()
181  detectorWcs = exposure.getWcs()
182  detectorCorners = detectorWcs.pixelToSky(geom.Box2D(detectorBBox).getCorners())
183  validPolygon = exposure.getInfo().getValidPolygon()
184  detectorPolygon = validPolygon if validPolygon else geom.Box2D(detectorBBox)
185 
186  availableCoaddRefs = dict()
187  overlappingArea = 0
188  for coaddRef in coaddExposureRefs:
189  dataId = coaddRef.dataId
190  patchWcs = skyMap[dataId['tract']].getWcs()
191  patchBBox = skyMap[dataId['tract']][dataId['patch']].getOuterBBox()
192  patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())
193  patchPolygon = afwGeom.Polygon(detectorWcs.skyToPixel(patchCorners))
194  if patchPolygon.intersection(detectorPolygon):
195  overlappingArea += patchPolygon.intersectionSingle(detectorPolygon).calculateArea()
196  if self.config.coaddName == 'dcr':
197  self.log.info("Using template input tract=%s, patch=%s, subfilter=%s",
198  dataId['tract'], dataId['patch'], dataId['subfilter'])
199  if dataId['patch'] in availableCoaddRefs:
200  availableCoaddRefs[dataId['patch']].append(coaddRef)
201  else:
202  availableCoaddRefs[dataId['patch']] = [coaddRef, ]
203  else:
204  self.log.info("Using template input tract=%s, patch=%s",
205  dataId['tract'], dataId['patch'])
206  availableCoaddRefs[dataId['patch']] = coaddRef
207 
208  if overlappingArea == 0:
209  templateExposure = None
210  pixGood = 0
211  self.log.warning("No overlapping template patches found")
212  else:
213  patchList = [tractInfo[patch] for patch in availableCoaddRefs.keys()]
214  templateExposure = self.runrun(tractInfo, patchList, detectorCorners, availableCoaddRefs,
215  visitInfo=exposure.getInfo().getVisitInfo())
216  # Count the number of pixels with the NO_DATA mask bit set
217  # counting NaN pixels is insufficient because pixels without data are often intepolated over)
218  pixNoData = np.count_nonzero(templateExposure.mask.array
219  & templateExposure.mask.getPlaneBitMask('NO_DATA'))
220  pixGood = templateExposure.getBBox().getArea() - pixNoData
221  self.log.info("template has %d good pixels (%.1f%%)", pixGood,
222  100*pixGood/templateExposure.getBBox().getArea())
223  return pipeBase.Struct(exposure=templateExposure, sources=None, area=pixGood)
224 
225  def getOverlapPatchList(self, exposure, skyMap):
226  """Select the relevant tract and its patches that overlap with the science exposure.
227 
228  Parameters
229  ----------
230  exposure : `lsst.afw.image.Exposure`
231  The science exposure to define the sky region of the template coadd.
232 
233  skyMap : `lsst.skymap.BaseSkyMap`
234  SkyMap object that corresponds to the template coadd.
235 
236  Returns
237  -------
238  result : `tuple` of
239  - ``tractInfo`` : `lsst.skymap.TractInfo`
240  The selected tract.
241  - ``patchList`` : `list` of `lsst.skymap.PatchInfo`
242  List of all overlap patches of the selected tract.
243  - ``skyCorners`` : `list` of `lsst.geom.SpherePoint`
244  Corners of the exposure in the sky in the order given by `lsst.geom.Box2D.getCorners`.
245  """
246  expWcs = exposure.getWcs()
247  expBoxD = geom.Box2D(exposure.getBBox())
248  expBoxD.grow(self.config.templateBorderSize)
249  ctrSkyPos = expWcs.pixelToSky(expBoxD.getCenter())
250  tractInfo = skyMap.findTract(ctrSkyPos)
251  self.log.info("Using skyMap tract %s", tractInfo.getId())
252  skyCorners = [expWcs.pixelToSky(pixPos) for pixPos in expBoxD.getCorners()]
253  patchList = tractInfo.findPatchList(skyCorners)
254 
255  if not patchList:
256  raise RuntimeError("No suitable tract found")
257 
258  self.log.info("Assembling %d coadd patches", len(patchList))
259  self.log.info("exposure dimensions=%s", exposure.getDimensions())
260 
261  return (tractInfo, patchList, skyCorners)
262 
263  def run(self, tractInfo, patchList, skyCorners, availableCoaddRefs,
264  sensorRef=None, visitInfo=None):
265  """Gen2 and gen3 shared code: determination of exposure dimensions and
266  copying of pixels from overlapping patch regions.
267 
268  Parameters
269  ----------
270  skyMap : `lsst.skymap.BaseSkyMap`
271  SkyMap object that corresponds to the template coadd.
272  tractInfo : `lsst.skymap.TractInfo`
273  The selected tract.
274  patchList : iterable of `lsst.skymap.patchInfo.PatchInfo`
275  Patches to consider for making the template exposure.
276  skyCorners : list of `lsst.geom.SpherePoint`
277  Sky corner coordinates to be covered by the template exposure.
278  availableCoaddRefs : `dict` [`int`]
279  Dictionary of spatially relevant retrieved coadd patches,
280  indexed by their sequential patch number. In Gen3 mode, values are
281  `lsst.daf.butler.DeferredDatasetHandle` and ``.get()`` is called,
282  in Gen2 mode, ``sensorRef.get(**coaddef)`` is called to retrieve the coadd.
283  sensorRef : `lsst.daf.persistence.ButlerDataRef`, Gen2 only
284  Butler data reference to get coadd data.
285  Must be `None` for Gen3.
286  visitInfo : `lsst.afw.image.VisitInfo`, Gen2 only
287  VisitInfo to make dcr model.
288 
289  Returns
290  -------
291  templateExposure: `lsst.afw.image.ExposureF`
292  The created template exposure.
293  """
294  coaddWcs = tractInfo.getWcs()
295 
296  # compute coadd bbox
297  coaddBBox = geom.Box2D()
298  for skyPos in skyCorners:
299  coaddBBox.include(coaddWcs.skyToPixel(skyPos))
300  coaddBBox = geom.Box2I(coaddBBox)
301  self.log.info("coadd dimensions=%s", coaddBBox.getDimensions())
302 
303  coaddExposure = afwImage.ExposureF(coaddBBox, coaddWcs)
304  coaddExposure.maskedImage.set(np.nan, afwImage.Mask.getPlaneBitMask("NO_DATA"), np.nan)
305  nPatchesFound = 0
306  coaddFilterLabel = None
307  coaddPsf = None
308  coaddPhotoCalib = None
309  for patchInfo in patchList:
310  patchNumber = tractInfo.getSequentialPatchIndex(patchInfo)
311  patchSubBBox = patchInfo.getOuterBBox()
312  patchSubBBox.clip(coaddBBox)
313  if patchNumber not in availableCoaddRefs:
314  self.log.warning("skip patch=%d; patch does not exist for this coadd", patchNumber)
315  continue
316  if patchSubBBox.isEmpty():
317  if isinstance(availableCoaddRefs[patchNumber], DeferredDatasetHandle):
318  tract = availableCoaddRefs[patchNumber].dataId['tract']
319  else:
320  tract = availableCoaddRefs[patchNumber]['tract']
321  self.log.info("skip tract=%d patch=%d; no overlapping pixels", tract, patchNumber)
322  continue
323 
324  if self.config.coaddName == 'dcr':
325  patchInnerBBox = patchInfo.getInnerBBox()
326  patchInnerBBox.clip(coaddBBox)
327  if np.min(patchInnerBBox.getDimensions()) <= 2*self.config.templateBorderSize:
328  self.log.info("skip tract=%(tract)s, patch=%(patch)s; too few pixels.",
329  availableCoaddRefs[patchNumber])
330  continue
331  self.log.info("Constructing DCR-matched template for patch %s",
332  availableCoaddRefs[patchNumber])
333 
334  if sensorRef:
335  dcrModel = DcrModel.fromDataRef(sensorRef,
336  self.config.effectiveWavelength,
337  self.config.bandwidth,
338  **availableCoaddRefs[patchNumber])
339  else:
340  dcrModel = DcrModel.fromQuantum(availableCoaddRefs[patchNumber],
341  self.config.effectiveWavelength,
342  self.config.bandwidth)
343  # The edge pixels of the DcrCoadd may contain artifacts due to missing data.
344  # Each patch has significant overlap, and the contaminated edge pixels in
345  # a new patch will overwrite good pixels in the overlap region from
346  # previous patches.
347  # Shrink the BBox to remove the contaminated pixels,
348  # but make sure it is only the overlap region that is reduced.
349  dcrBBox = geom.Box2I(patchSubBBox)
350  dcrBBox.grow(-self.config.templateBorderSize)
351  dcrBBox.include(patchInnerBBox)
352  coaddPatch = dcrModel.buildMatchedExposure(bbox=dcrBBox,
353  wcs=coaddWcs,
354  visitInfo=visitInfo)
355  else:
356  if sensorRef is None:
357  # Gen3
358  coaddPatch = availableCoaddRefs[patchNumber].get()
359  else:
360  # Gen2
361  coaddPatch = sensorRef.get(**availableCoaddRefs[patchNumber])
362  nPatchesFound += 1
363 
364  # Gen2 get() seems to clip based on bbox kwarg but we removed bbox
365  # calculation from caller code. Gen3 also does not do this.
366  overlapBox = coaddPatch.getBBox()
367  overlapBox.clip(coaddBBox)
368  coaddExposure.maskedImage.assign(coaddPatch.maskedImage[overlapBox], overlapBox)
369 
370  if coaddFilterLabel is None:
371  coaddFilterLabel = coaddPatch.getFilterLabel()
372 
373  # Retrieve the PSF for this coadd tract, if not already retrieved
374  if coaddPsf is None and coaddPatch.hasPsf():
375  coaddPsf = coaddPatch.getPsf()
376 
377  # Retrieve the calibration for this coadd tract, if not already retrieved
378  if coaddPhotoCalib is None:
379  coaddPhotoCalib = coaddPatch.getPhotoCalib()
380 
381  if coaddPhotoCalib is None:
382  raise RuntimeError("No coadd PhotoCalib found!")
383  if nPatchesFound == 0:
384  raise RuntimeError("No patches found!")
385  if coaddPsf is None:
386  raise RuntimeError("No coadd Psf found!")
387 
388  coaddExposure.setPhotoCalib(coaddPhotoCalib)
389  coaddExposure.setPsf(coaddPsf)
390  coaddExposure.setFilterLabel(coaddFilterLabel)
391  return coaddExposure
392 
394  """Return coadd name for given task config
395 
396  Returns
397  -------
398  CoaddDatasetName : `string`
399 
400  TODO: This nearly duplicates a method in CoaddBaseTask (DM-11985)
401  """
402  warpType = self.config.warpType
403  suffix = "" if warpType == "direct" else warpType[0].upper() + warpType[1:]
404  return self.config.coaddName + "Coadd" + suffix
405 
406 
407 class GetCalexpAsTemplateConfig(pexConfig.Config):
408  doAddCalexpBackground = pexConfig.Field(
409  dtype=bool,
410  default=True,
411  doc="Add background to calexp before processing it."
412  )
413 
414 
415 class GetCalexpAsTemplateTask(pipeBase.Task):
416  """Subtask to retrieve calexp of the same ccd number as the science image SensorRef
417  for use as an image difference template. Only gen2 supported.
418 
419  To be run as a subtask by pipe.tasks.ImageDifferenceTask.
420  Intended for use with simulations and surveys that repeatedly visit the same pointing.
421  This code was originally part of Winter2013ImageDifferenceTask.
422  """
423 
424  ConfigClass = GetCalexpAsTemplateConfig
425  _DefaultName = "GetCalexpAsTemplateTask"
426 
427  def run(self, exposure, sensorRef, templateIdList):
428  """Return a calexp exposure with based on input sensorRef.
429 
430  Construct a dataId based on the sensorRef.dataId combined
431  with the specifications from the first dataId in templateIdList
432 
433  Parameters
434  ----------
435  exposure : `lsst.afw.image.Exposure`
436  exposure (unused)
437  sensorRef : `list` of `lsst.daf.persistence.ButlerDataRef`
438  Data reference of the calexp(s) to subtract from.
439  templateIdList : `list` of `lsst.daf.persistence.ButlerDataRef`
440  Data reference of the template calexp to be subtraced.
441  Can be incomplete, fields are initialized from `sensorRef`.
442  If there are multiple items, only the first one is used.
443 
444  Returns
445  -------
446  result : `struct`
447 
448  return a pipeBase.Struct:
449 
450  - ``exposure`` : a template calexp
451  - ``sources`` : source catalog measured on the template
452  """
453 
454  if len(templateIdList) == 0:
455  raise RuntimeError("No template data reference supplied.")
456  if len(templateIdList) > 1:
457  self.log.warning("Multiple template data references supplied. Using the first one only.")
458 
459  templateId = sensorRef.dataId.copy()
460  templateId.update(templateIdList[0])
461 
462  self.log.info("Fetching calexp (%s) as template.", templateId)
463 
464  butler = sensorRef.getButler()
465  template = butler.get(datasetType="calexp", dataId=templateId)
466  if self.config.doAddCalexpBackground:
467  templateBg = butler.get(datasetType="calexpBackground", dataId=templateId)
468  mi = template.getMaskedImage()
469  mi += templateBg.getImage()
470 
471  if not template.hasPsf():
472  raise pipeBase.TaskError("Template has no psf")
473 
474  templateSources = butler.get(datasetType="src", dataId=templateId)
475  return pipeBase.Struct(exposure=template,
476  sources=templateSources)
477 
478  def runDataRef(self, *args, **kwargs):
479  return self.runrun(*args, **kwargs)
480 
481  def runQuantum(self, **kwargs):
482  raise NotImplementedError("Calexp template is not supported with gen3 middleware")
483 
484 
485 class GetMultiTractCoaddTemplateConnections(pipeBase.PipelineTaskConnections,
486  dimensions=("instrument", "visit", "detector", "skymap"),
487  defaultTemplates={"coaddName": "goodSeeing",
488  "warpTypeSuffix": "",
489  "fakesType": ""}):
490  bbox = pipeBase.connectionTypes.Input(
491  doc="BBoxes of calexp used determine geometry of output template",
492  name="{fakesType}calexp.bbox",
493  storageClass="Box2I",
494  dimensions=("instrument", "visit", "detector"),
495  )
496  wcs = pipeBase.connectionTypes.Input(
497  doc="WCSs of calexps that we want to fetch the template for",
498  name="{fakesType}calexp.wcs",
499  storageClass="Wcs",
500  dimensions=("instrument", "visit", "detector"),
501  )
502  skyMap = pipeBase.connectionTypes.Input(
503  doc="Input definition of geometry/bbox and projection/wcs for template exposures",
504  name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
505  dimensions=("skymap", ),
506  storageClass="SkyMap",
507  )
508  # TODO DM-31292: Add option to use global external wcs from jointcal
509  # Needed for DRP HSC
510  coaddExposures = pipeBase.connectionTypes.Input(
511  doc="Input template to match and subtract from the exposure",
512  dimensions=("tract", "patch", "skymap", "band"),
513  storageClass="ExposureF",
514  name="{fakesType}{coaddName}Coadd{warpTypeSuffix}",
515  multiple=True,
516  deferLoad=True
517  )
518  outputExposure = pipeBase.connectionTypes.Output(
519  doc="Warped template used to create `subtractedExposure`.",
520  dimensions=("instrument", "visit", "detector"),
521  storageClass="ExposureF",
522  name="{fakesType}{coaddName}Diff_templateExp{warpTypeSuffix}",
523  )
524 
525 
526 class GetMultiTractCoaddTemplateConfig(pipeBase.PipelineTaskConfig, GetCoaddAsTemplateConfig,
527  pipelineConnections=GetMultiTractCoaddTemplateConnections):
528  warp = pexConfig.ConfigField(
529  dtype=afwMath.Warper.ConfigClass,
530  doc="warper configuration",
531  )
532  coaddPsf = pexConfig.ConfigField(
533  doc="Configuration for CoaddPsf",
534  dtype=CoaddPsfConfig,
535  )
536 
537  def setDefaults(self):
538  self.warp.warpingKernelName = 'lanczos5'
539  self.coaddPsf.warpingKernelName = 'lanczos5'
540 
541 
542 class GetMultiTractCoaddTemplateTask(pipeBase.PipelineTask):
543  ConfigClass = GetMultiTractCoaddTemplateConfig
544  _DefaultName = "getMultiTractCoaddTemplateTask"
545 
546  def __init__(self, *args, **kwargs):
547  super().__init__(*args, **kwargs)
548  self.warper = afwMath.Warper.fromConfig(self.config.warp)
549 
550  def runQuantum(self, butlerQC, inputRefs, outputRefs):
551  # Read in all inputs.
552  inputs = butlerQC.get(inputRefs)
553  inputs['coaddExposures'] = self.getOverlappingExposures(inputs)
554  # SkyMap only needed for filtering without
555  inputs.pop('skyMap')
556  outputs = self.run(**inputs)
557  butlerQC.put(outputs, outputRefs)
558 
559  def getOverlappingExposures(self, inputs):
560  """Return list of coaddExposure DeferredDatasetHandles that overlap detector
561 
562  The spatial index in the registry has generous padding and often supplies
563  patches near, but not directly overlapping the detector.
564  Filters inputs so that we don't have to read in all input coadds.
565 
566  Parameters
567  ----------
568  inputs : `dict` of task Inputs
569 
570  Returns
571  -------
572  coaddExposures : list of elements of type
573  `lsst.daf.butler.DeferredDatasetHandle` of
574  `lsst.afw.image.Exposure`
575 
576  Raises
577  ------
578  NoWorkFound
579  Raised if no patches overlap the input detector bbox
580  """
581  # Check that the patches actually overlap the detector
582  # Exposure's validPolygon would be more accurate
583  detectorPolygon = geom.Box2D(inputs['bbox'])
584  overlappingArea = 0
585  coaddExposureList = []
586  for coaddRef in inputs['coaddExposures']:
587  dataId = coaddRef.dataId
588  patchWcs = inputs['skyMap'][dataId['tract']].getWcs()
589  patchBBox = inputs['skyMap'][dataId['tract']][dataId['patch']].getOuterBBox()
590  patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())
591  patchPolygon = afwGeom.Polygon(inputs['wcs'].skyToPixel(patchCorners))
592  if patchPolygon.intersection(detectorPolygon):
593  overlappingArea += patchPolygon.intersectionSingle(detectorPolygon).calculateArea()
594  self.log.info("Using template input tract=%s, patch=%s" %
595  (dataId['tract'], dataId['patch']))
596  coaddExposureList.append(coaddRef)
597 
598  if not overlappingArea:
599  raise pipeBase.NoWorkFound('No patches overlap detector')
600 
601  return coaddExposureList
602 
603  def run(self, coaddExposures, bbox, wcs):
604  """Warp coadds from multiple tracts to form a template for image diff.
605 
606  Where the tracts overlap, the resulting template image is averaged.
607  The PSF on the template is created by combining the CoaddPsf on each
608  template image into a meta-CoaddPsf.
609 
610  Parameters
611  ----------
612  coaddExposures: list of DeferredDatasetHandle to `lsst.afw.image.Exposure`
613  Coadds to be mosaicked
614  bbox : `lsst.geom.Box2I`
615  Template Bounding box of the detector geometry onto which to
616  resample the coaddExposures
617  wcs : `lsst.afw.geom.SkyWcs`
618  Template WCS onto which to resample the coaddExposures
619 
620  Returns
621  -------
622  result : `struct`
623  return a pipeBase.Struct:
624  - ``outputExposure`` : a template coadd exposure assembled out of patches
625 
626 
627  Raises
628  ------
629  NoWorkFound
630  Raised if no patches overlatp the input detector bbox
631 
632  """
633  # Table for CoaddPSF
634  tractsSchema = afwTable.ExposureTable.makeMinimalSchema()
635  tractKey = tractsSchema.addField('tract', type=np.int32, doc='Which tract')
636  patchKey = tractsSchema.addField('patch', type=np.int32, doc='Which patch')
637  weightKey = tractsSchema.addField('weight', type=float, doc='Weight for each tract, should be 1')
638  tractsCatalog = afwTable.ExposureCatalog(tractsSchema)
639 
640  finalWcs = wcs
641  bbox.grow(self.config.templateBorderSize)
642  finalBBox = bbox
643 
644  nPatchesFound = 0
645  maskedImageList = []
646  weightList = []
647 
648  for coaddExposure in coaddExposures:
649  coaddPatch = coaddExposure.get()
650 
651  # warp to detector WCS
652  xyTransform = afwGeom.makeWcsPairTransform(coaddPatch.getWcs(), finalWcs)
653  psfWarped = WarpedPsf(coaddPatch.getPsf(), xyTransform)
654  warped = self.warper.warpExposure(finalWcs, coaddPatch, maxBBox=finalBBox)
655 
656  # Check if warped image is viable
657  if not np.any(np.isfinite(warped.image.array)):
658  self.log.info("No overlap for warped %s. Skipping" % coaddExposure.ref.dataId)
659  continue
660 
661  warped.setPsf(psfWarped)
662 
663  exp = afwImage.ExposureF(finalBBox, finalWcs)
664  exp.maskedImage.set(np.nan, afwImage.Mask.getPlaneBitMask("NO_DATA"), np.nan)
665  exp.maskedImage.assign(warped.maskedImage, warped.getBBox())
666 
667  maskedImageList.append(exp.maskedImage)
668  weightList.append(1)
669  record = tractsCatalog.addNew()
670  record.setPsf(psfWarped)
671  record.setWcs(finalWcs)
672  record.setPhotoCalib(coaddPatch.getPhotoCalib())
673  record.setBBox(warped.getBBox())
674  record.set(tractKey, coaddExposure.ref.dataId['tract'])
675  record.set(patchKey, coaddExposure.ref.dataId['patch'])
676  record.set(weightKey, 1.)
677  nPatchesFound += 1
678 
679  if nPatchesFound == 0:
680  raise pipeBase.NoWorkFound("No patches found to overlap detector")
681 
682  # Combine images from individual patches together
683  statsFlags = afwMath.stringToStatisticsProperty('MEAN')
684  statsCtrl = afwMath.StatisticsControl()
685  statsCtrl.setNanSafe(True)
686  statsCtrl.setWeighted(True)
687  statsCtrl.setCalcErrorFromInputVariance(True)
688 
689  templateExposure = afwImage.ExposureF(finalBBox, finalWcs)
690  templateExposure.maskedImage.set(np.nan, afwImage.Mask.getPlaneBitMask("NO_DATA"), np.nan)
691  xy0 = templateExposure.getXY0()
692  # Do not mask any values
693  templateExposure.maskedImage = afwMath.statisticsStack(maskedImageList, statsFlags, statsCtrl,
694  weightList, clipped=0, maskMap=[])
695  templateExposure.maskedImage.setXY0(xy0)
696 
697  # CoaddPsf centroid not only must overlap image, but must overlap the part of
698  # image with data. Use centroid of region with data
699  boolmask = templateExposure.mask.array & templateExposure.mask.getPlaneBitMask('NO_DATA') == 0
700  maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel))
701  centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
702 
703  ctrl = self.config.coaddPsf.makeControl()
704  coaddPsf = CoaddPsf(tractsCatalog, finalWcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize)
705  if coaddPsf is None:
706  raise RuntimeError("CoaddPsf could not be constructed")
707 
708  templateExposure.setPsf(coaddPsf)
709  templateExposure.setFilterLabel(coaddPatch.getFilterLabel())
710  templateExposure.setPhotoCalib(coaddPatch.getPhotoCalib())
711  return pipeBase.Struct(outputExposure=templateExposure)
Cartesian polygons.
Definition: Polygon.h:59
Pass parameters to a Statistics object.
Definition: Statistics.h:92
Custom catalog class for ExposureRecord/Table.
Definition: Exposure.h:311
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
An integer coordinate rectangle.
Definition: Box.h:55
def run(self, exposure, sensorRef, templateIdList)
Definition: getTemplate.py:427
def runDataRef(self, exposure, sensorRef, templateIdList=None)
Definition: getTemplate.py:104
def run(self, tractInfo, patchList, skyCorners, availableCoaddRefs, sensorRef=None, visitInfo=None)
Definition: getTemplate.py:264
def runQuantum(self, exposure, butlerQC, skyMapRef, coaddExposureRefs)
Definition: getTemplate.py:149
CoaddPsf is the Psf derived to be used for non-PSF-matched Coadd images.
Definition: CoaddPsf.h:58
A Psf class that maps an arbitrary Psf through a coordinate transformation.
Definition: WarpedPsf.h:52
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
Definition: SkyWcs.cc:146
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
Definition: Statistics.cc:738
std::shared_ptr< lsst::afw::image::Image< PixelT > > statisticsStack(std::vector< std::shared_ptr< lsst::afw::image::Image< PixelT >>> &images, Property flags, StatisticsControl const &sctrl=StatisticsControl(), std::vector< lsst::afw::image::VariancePixel > const &wvector=std::vector< lsst::afw::image::VariancePixel >(0))
A function to compute some statistics of a stack of Images.
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.
def run(self, coaddExposures, bbox, wcs)
Definition: getTemplate.py:603
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.