Loading [MathJax]/extensions/tex2jax.js
LSST Applications g0fba68d861+6234375589,g1e78f5e6d3+628c932f15,g1fd858c14a+eb8d917efb,g35bb328faa+fcb1d3bbc8,g4af146b050+4faa9dad44,g4d2262a081+3c0ac0cdae,g4e0f332c67+8616b824a5,g53246c7159+fcb1d3bbc8,g5a012ec0e7+d65fd7031a,g60b5630c4e+042d43a120,g67b6fd64d1+c0248a1c13,g6ea0df0560+042d43a120,g78460c75b0+2f9a1b4bcd,g786e29fd12+cf7ec2a62a,g7b71ed6315+fcb1d3bbc8,g87b7deb4dc+7fa294d175,g8852436030+40f6ec51d1,g89139ef638+c0248a1c13,g9125e01d80+fcb1d3bbc8,g94187f82dc+042d43a120,g989de1cb63+c0248a1c13,g9f33ca652e+da0da5ecef,g9f7030ddb1+682810b470,ga2b97cdc51+042d43a120,gabe3b4be73+1e0a283bba,gabf8522325+83c19109ce,gb1101e3267+5921e058d2,gb58c049af0+f03b321e39,gb89ab40317+c0248a1c13,gcf25f946ba+40f6ec51d1,gd6cbbdb0b4+d9e8db455e,gd9a9a58781+fcb1d3bbc8,gdabf7e867e+a3799d3da4,gde0f65d7ad+4a3db72839,ge278dab8ac+4ce6343b44,ge410e46f29+c0248a1c13,gf67bdafdda+c0248a1c13,gfe06eef73a+95f9f0e40c,v29.0.0.rc3
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
lsst.ip.diffim.getTemplate Namespace Reference

Classes

class  GetDcrTemplateConnections
 
class  GetTemplateConnections
 

Functions

 run (self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter)
 
 _makeExposureCatalog (self, exposureRefs, dataIds)
 
 _merge (self, maskedImages, bbox, wcs)
 
 _makePsf (self, template, catalog, wcs)
 
 checkPatchList (self, patchList)
 
 getDcrModel (self, patchList, coaddRefs, visitInfo)
 
 _selectDataRef (coaddRef, tract, patch)
 

Variables

 detectorPolygon = geom.Box2D(bbox)
 
int overlappingArea = 0
 
 coaddExposures = collections.defaultdict(list)
 
 dataIds = collections.defaultdict(list)
 
 dataId = coaddRef.dataId
 
 patchWcs = skymap[dataId['tract']].getWcs()
 
 patchBBox = skymap[dataId['tract']][dataId['patch']].getOuterBBox()
 
 patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())
 
 patchPolygon = afwGeom.Polygon(wcs.skyToPixel(patchCorners))
 
 patchList = dict()
 

Function Documentation

◆ _makeExposureCatalog()

lsst.ip.diffim.getTemplate._makeExposureCatalog ( self,
exposureRefs,
dataIds )
protected
Make an exposure catalog for one tract.

Parameters
----------
exposureRefs : `list` of [`lsst.daf.butler.DeferredDatasetHandle` of \
                `lsst.afw.image.Exposure`]
    Exposures to include in the catalog.
dataIds : `list` [`lsst.daf.butler.DataCoordinate`]
    Data ids of each of the included exposures; must have "tract" and
    "patch" entries.

Returns
-------
images : `dict` [`lsst.afw.image.MaskedImage`]
    MaskedImages of each of the input exposures, for warping.
catalog : `lsst.afw.table.ExposureCatalog`
    Catalog of metadata for each exposure
totalBox : `lsst.geom.Box2I`
    The union of the bounding boxes of all the input exposures.

Definition at line 366 of file getTemplate.py.

366 def _makeExposureCatalog(self, exposureRefs, dataIds):
367 """Make an exposure catalog for one tract.
368
369 Parameters
370 ----------
371 exposureRefs : `list` of [`lsst.daf.butler.DeferredDatasetHandle` of \
372 `lsst.afw.image.Exposure`]
373 Exposures to include in the catalog.
374 dataIds : `list` [`lsst.daf.butler.DataCoordinate`]
375 Data ids of each of the included exposures; must have "tract" and
376 "patch" entries.
377
378 Returns
379 -------
380 images : `dict` [`lsst.afw.image.MaskedImage`]
381 MaskedImages of each of the input exposures, for warping.
382 catalog : `lsst.afw.table.ExposureCatalog`
383 Catalog of metadata for each exposure
384 totalBox : `lsst.geom.Box2I`
385 The union of the bounding boxes of all the input exposures.
386 """
387 catalog = afwTable.ExposureCatalog(self.schema)
388 catalog.reserve(len(exposureRefs))
389 exposures = (exposureRef.get() for exposureRef in exposureRefs)
390 images = {}
391 totalBox = geom.Box2I()
392
393 for coadd, dataId in zip(exposures, dataIds):
394 images[dataId] = coadd.maskedImage
395 bbox = coadd.getBBox()
396 totalBox = totalBox.expandedTo(bbox)
397 record = catalog.addNew()
398 record.setPsf(coadd.psf)
399 record.setWcs(coadd.wcs)
400 record.setPhotoCalib(coadd.photoCalib)
401 record.setBBox(bbox)
402 record.setValidPolygon(afwGeom.Polygon(geom.Box2D(bbox).getCorners()))
403 record.set("tract", dataId["tract"])
404 record.set("patch", dataId["patch"])
405 # Weight is used by CoaddPsf, but the PSFs from overlapping patches
406 # should be very similar, so this value mostly shouldn't matter.
407 record.set("weight", 1)
408
409 return images, catalog, totalBox
410
Cartesian polygons.
Definition Polygon.h:59
A floating-point coordinate rectangle geometry.
Definition Box.h:413
An integer coordinate rectangle.
Definition Box.h:55

◆ _makePsf()

lsst.ip.diffim.getTemplate._makePsf ( self,
template,
catalog,
wcs )
protected
Return a PSF containing the PSF at each of the input regions.

Note that although this includes all the exposures from the catalog,
the PSF knows which part of the template the inputs came from, so when
evaluated at a given position it will not include inputs that never
went in to those pixels.

Parameters
----------
template : `lsst.afw.image.Exposure`
    Generated template the PSF is for.
catalog : `lsst.afw.table.ExposureCatalog`
    Catalog of exposures that went into the template that contains all
    of the input PSFs.
wcs : `lsst.afw.geom.SkyWcs`
    WCS of the template, to warp the PSFs to.

Returns
-------
coaddPsf : `lsst.meas.algorithms.CoaddPsf`
    The meta-psf constructed from all of the input catalogs.

Definition at line 483 of file getTemplate.py.

483 def _makePsf(self, template, catalog, wcs):
484 """Return a PSF containing the PSF at each of the input regions.
485
486 Note that although this includes all the exposures from the catalog,
487 the PSF knows which part of the template the inputs came from, so when
488 evaluated at a given position it will not include inputs that never
489 went in to those pixels.
490
491 Parameters
492 ----------
493 template : `lsst.afw.image.Exposure`
494 Generated template the PSF is for.
495 catalog : `lsst.afw.table.ExposureCatalog`
496 Catalog of exposures that went into the template that contains all
497 of the input PSFs.
498 wcs : `lsst.afw.geom.SkyWcs`
499 WCS of the template, to warp the PSFs to.
500
501 Returns
502 -------
503 coaddPsf : `lsst.meas.algorithms.CoaddPsf`
504 The meta-psf constructed from all of the input catalogs.
505 """
506 # CoaddPsf centroid not only must overlap image, but must overlap the
507 # part of image with data. Use centroid of region with data.
508 boolmask = template.mask.array & template.mask.getPlaneBitMask('NO_DATA') == 0
509 maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel))
510 centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
511
512 ctrl = self.config.coaddPsf.makeControl()
513 coaddPsf = CoaddPsf(catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize)
514 return coaddPsf
515
516

◆ _merge()

lsst.ip.diffim.getTemplate._merge ( self,
maskedImages,
bbox,
wcs )
protected
Merge the images that came from one tract into one larger image,
ignoring NaN pixels and non-finite variance pixels from individual
exposures.

Parameters
----------
maskedImages : `dict` [`lsst.afw.image.MaskedImage` or
                       `lsst.afw.image.Exposure`]
    Images to be merged into one larger bounding box.
bbox : `lsst.geom.Box2I`
    Bounding box defining the image to merge into.
wcs : `lsst.afw.geom.SkyWcs`
    WCS of all of the input images to set on the output image.

Returns
-------
merged : `lsst.afw.image.MaskedImage`
    Merged image with all of the inputs at their respective bbox
    positions.
count : `int`
    Count of the number of good pixels (those with positive weights)
    in the merged image.
included : `list` [`int`]
    List of indexes of patches that were included in the merged
    result, to be used to trim the exposure catalog.

Definition at line 411 of file getTemplate.py.

411 def _merge(self, maskedImages, bbox, wcs):
412 """Merge the images that came from one tract into one larger image,
413 ignoring NaN pixels and non-finite variance pixels from individual
414 exposures.
415
416 Parameters
417 ----------
418 maskedImages : `dict` [`lsst.afw.image.MaskedImage` or
419 `lsst.afw.image.Exposure`]
420 Images to be merged into one larger bounding box.
421 bbox : `lsst.geom.Box2I`
422 Bounding box defining the image to merge into.
423 wcs : `lsst.afw.geom.SkyWcs`
424 WCS of all of the input images to set on the output image.
425
426 Returns
427 -------
428 merged : `lsst.afw.image.MaskedImage`
429 Merged image with all of the inputs at their respective bbox
430 positions.
431 count : `int`
432 Count of the number of good pixels (those with positive weights)
433 in the merged image.
434 included : `list` [`int`]
435 List of indexes of patches that were included in the merged
436 result, to be used to trim the exposure catalog.
437 """
438 merged = afwImage.ExposureF(bbox, wcs)
439 weights = afwImage.ImageF(bbox)
440 included = [] # which patches were included in the result
441 for i, (dataId, maskedImage) in enumerate(maskedImages.items()):
442 # Only merge into the trimmed box, to save memory
443 clippedBox = geom.Box2I(maskedImage.getBBox())
444 clippedBox.clip(bbox)
445 if clippedBox.area == 0:
446 self.log.debug("%s does not overlap template region.", dataId)
447 continue # nothing in this image overlaps the output
448 maskedImage = maskedImage.subset(clippedBox)
449 # Catch both zero-value and NaN variance plane pixels
450 good = (maskedImage.variance.array > 0) & (np.isfinite(maskedImage.variance.array))
451 weight = maskedImage.variance.array[good]**(-0.5)
452 bad = np.isnan(maskedImage.image.array) | ~good
453 # Note that modifying the patch MaskedImage in place is fine;
454 # we're throwing it away at the end anyway.
455 maskedImage.image.array[bad] = 0.0
456 maskedImage.variance.array[bad] = 0.0
457 # Reset mask, too, since these pixels don't contribute to sum.
458 maskedImage.mask.array[bad] = 0
459 # Cannot use `merged.maskedImage *= weight` because that operator
460 # multiplies the variance by the weight twice; in this case
461 # `weight` are the exact values we want to scale by.
462 maskedImage.image.array[good] *= weight
463 maskedImage.variance.array[good] *= weight
464 weights[clippedBox].array[good] += weight
465 # Free memory before creating new large arrays
466 del weight
467 merged.maskedImage[clippedBox] += maskedImage
468 included.append(i)
469
470 good = weights.array > 0
471
472 # Cannot use `merged.maskedImage /= weights` because that
473 # operator divides the variance by the weight twice; in this case
474 # `weights` are the exact values we want to scale by.
475 weights = weights.array[good]
476 merged.image.array[good] /= weights
477 merged.variance.array[good] /= weights
478
479 merged.mask.array[~good] |= merged.mask.getPlaneBitMask("NO_DATA")
480
481 return merged, good.sum(), included
482

◆ _selectDataRef()

lsst.ip.diffim.getTemplate._selectDataRef ( coaddRef,
tract,
patch )
protected

Definition at line 732 of file getTemplate.py.

732def _selectDataRef(coaddRef, tract, patch):
733 condition = (coaddRef.dataId['tract'] == tract) & (coaddRef.dataId['patch'] == patch)
734 return condition

◆ checkPatchList()

lsst.ip.diffim.getTemplate.checkPatchList ( self,
patchList )
Check that all of the DcrModel subfilters are present for each
patch.

Parameters
----------
patchList : `dict`
    Dict of the patches containing valid data for each tract.

Raises
------
RuntimeError
    If the number of exposures found for a patch does not match the
    number of subfilters.

Definition at line 679 of file getTemplate.py.

679 def checkPatchList(self, patchList):
680 """Check that all of the DcrModel subfilters are present for each
681 patch.
682
683 Parameters
684 ----------
685 patchList : `dict`
686 Dict of the patches containing valid data for each tract.
687
688 Raises
689 ------
690 RuntimeError
691 If the number of exposures found for a patch does not match the
692 number of subfilters.
693 """
694 for tract in patchList:
695 for patch in set(patchList[tract]):
696 if patchList[tract].count(patch) != self.config.numSubfilters:
697 raise RuntimeError("Invalid number of DcrModel subfilters found: %d vs %d expected",
698 patchList[tract].count(patch), self.config.numSubfilters)
699

◆ getDcrModel()

lsst.ip.diffim.getTemplate.getDcrModel ( self,
patchList,
coaddRefs,
visitInfo )
Build DCR-matched coadds from a list of exposure references.

Parameters
----------
patchList : `dict`
    Dict of the patches containing valid data for each tract.
coaddRefs : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
    Data references to `~lsst.afw.image.Exposure` representing
    DcrModels that overlap the detector.
visitInfo : `lsst.afw.image.VisitInfo`
    Metadata for the science image.

Returns
-------
coaddExposures : `list` [`lsst.afw.image.Exposure`]
    Coadd exposures that overlap the detector.

Definition at line 700 of file getTemplate.py.

700 def getDcrModel(self, patchList, coaddRefs, visitInfo):
701 """Build DCR-matched coadds from a list of exposure references.
702
703 Parameters
704 ----------
705 patchList : `dict`
706 Dict of the patches containing valid data for each tract.
707 coaddRefs : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
708 Data references to `~lsst.afw.image.Exposure` representing
709 DcrModels that overlap the detector.
710 visitInfo : `lsst.afw.image.VisitInfo`
711 Metadata for the science image.
712
713 Returns
714 -------
715 coaddExposures : `list` [`lsst.afw.image.Exposure`]
716 Coadd exposures that overlap the detector.
717 """
718 coaddExposures = collections.defaultdict(list)
719 for tract in patchList:
720 for patch in set(patchList[tract]):
721 coaddRefList = [coaddRef for coaddRef in coaddRefs
722 if _selectDataRef(coaddRef, tract, patch)]
723
724 dcrModel = DcrModel.fromQuantum(coaddRefList,
725 self.config.effectiveWavelength,
726 self.config.bandwidth,
727 self.config.numSubfilters)
728 coaddExposures[tract].append(dcrModel.buildMatchedExposure(visitInfo=visitInfo))
729 return coaddExposures
730
731

◆ run()

lsst.ip.diffim.getTemplate.run ( self,
* ,
coaddExposureHandles,
bbox,
wcs,
dataIds,
physical_filter )
Warp coadds from multiple tracts and patches to form a template to
subtract from a science image.

Tract and patch overlap regions are combined by a variance-weighted
average, and the variance planes are combined with the same weights,
not added in quadrature; the overlap regions are not statistically
independent, because they're derived from the same original data.
The PSF on the template is created by combining the CoaddPsf on each
template image into a meta-CoaddPsf.

Parameters
----------
coaddExposureHandles : `dict` [`int`,  `list` of \
                  [`lsst.daf.butler.DeferredDatasetHandle` of \
                   `lsst.afw.image.Exposure`]]
    Coadds to be mosaicked, indexed on tract id.
bbox : `lsst.geom.Box2I`
    Template Bounding box of the detector geometry onto which to
    resample the ``coaddExposureHandles``. Modified in-place to include the
    template border.
wcs : `lsst.afw.geom.SkyWcs`
    Template WCS onto which to resample the ``coaddExposureHandles``.
dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
    Record of the tract and patch of each coaddExposure, indexed on
    tract id.
physical_filter : `str`
    Physical filter of the science image.

Returns
-------
result : `lsst.pipe.base.Struct`
   A struct with attributes:

   ``template``
       A template coadd exposure assembled out of patches
       (`lsst.afw.image.ExposureF`).

Raises
------
NoWorkFound
    If no coadds are found with sufficient un-masked pixels.

Definition at line 229 of file getTemplate.py.

229 def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter):
230 """Warp coadds from multiple tracts and patches to form a template to
231 subtract from a science image.
232
233 Tract and patch overlap regions are combined by a variance-weighted
234 average, and the variance planes are combined with the same weights,
235 not added in quadrature; the overlap regions are not statistically
236 independent, because they're derived from the same original data.
237 The PSF on the template is created by combining the CoaddPsf on each
238 template image into a meta-CoaddPsf.
239
240 Parameters
241 ----------
242 coaddExposureHandles : `dict` [`int`, `list` of \
243 [`lsst.daf.butler.DeferredDatasetHandle` of \
244 `lsst.afw.image.Exposure`]]
245 Coadds to be mosaicked, indexed on tract id.
246 bbox : `lsst.geom.Box2I`
247 Template Bounding box of the detector geometry onto which to
248 resample the ``coaddExposureHandles``. Modified in-place to include the
249 template border.
250 wcs : `lsst.afw.geom.SkyWcs`
251 Template WCS onto which to resample the ``coaddExposureHandles``.
252 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
253 Record of the tract and patch of each coaddExposure, indexed on
254 tract id.
255 physical_filter : `str`
256 Physical filter of the science image.
257
258 Returns
259 -------
260 result : `lsst.pipe.base.Struct`
261 A struct with attributes:
262
263 ``template``
264 A template coadd exposure assembled out of patches
265 (`lsst.afw.image.ExposureF`).
266
267 Raises
268 ------
269 NoWorkFound
270 If no coadds are found with sufficient un-masked pixels.
271 """
272 band, photoCalib = self._checkInputs(dataIds, coaddExposureHandles)
273
274 bbox.grow(self.config.templateBorderSize)
275
276 warped = {}
277 catalogs = []
278 for tract in coaddExposureHandles:
279 maskedImages, catalog, totalBox = self._makeExposureCatalog(coaddExposureHandles[tract],
280 dataIds[tract])
281 warpedBox = computeWarpedBBox(catalog[0].wcs, bbox, wcs)
282 warpedBox.grow(5) # to ensure we catch all relevant input pixels
283 # Combine images from individual patches together.
284 unwarped, count, included = self._merge(maskedImages, warpedBox, catalog[0].wcs)
285 # Delete `maskedImages` after combining into one large image to reduce peak memory use
286 del maskedImages
287 if count == 0:
288 self.log.info("No valid pixels from coadd patches in tract %s; not including in output.",
289 tract)
290 continue
291 warpedBox.clip(totalBox)
292 potentialInput = self.warper.warpExposure(wcs, unwarped.subset(warpedBox), destBBox=bbox)
293
294 # Delete the single large `unwarped` image after warping to reduce peak memory use
295 del unwarped
296 if np.all(potentialInput.mask.array & potentialInput.mask.getPlaneBitMask("NO_DATA")):
297 self.log.info("No overlap from coadd patches in tract %s; not including in output.", tract)
298 continue
299
300 # Trim the exposure catalog to just the patches that were used.
301 tempCatalog = afwTable.ExposureCatalog(self.schema)
302 tempCatalog.reserve(len(included))
303 for i in included:
304 tempCatalog.append(catalog[i])
305 catalogs.append(tempCatalog)
306 warped[tract] = potentialInput.maskedImage
307
308 if len(warped) == 0:
309 raise pipeBase.NoWorkFound("No patches found to overlap science exposure.")
310 # At this point, all entries will be valid, so we can ignore included.
311 template, count, _ = self._merge(warped, bbox, wcs)
312 if count == 0:
313 raise pipeBase.NoWorkFound("No valid pixels in warped template.")
314
315 # Make a single catalog containing all the inputs that were accepted.
316 catalog = afwTable.ExposureCatalog(self.schema)
317 catalog.reserve(sum([len(c) for c in catalogs]))
318 for c in catalogs:
319 catalog.extend(c)
320
321 template.setPsf(self._makePsf(template, catalog, wcs))
322 template.setFilter(afwImage.FilterLabel(band, physical_filter))
323 template.setPhotoCalib(photoCalib)
324 return pipeBase.Struct(template=template)
325
A group of labels for a filter in an exposure or coadd.
Definition FilterLabel.h:58

Variable Documentation

◆ coaddExposures

lsst.ip.diffim.getTemplate.coaddExposures = collections.defaultdict(list)

Definition at line 207 of file getTemplate.py.

◆ dataId

lsst.ip.diffim.getTemplate.dataId = coaddRef.dataId

Definition at line 211 of file getTemplate.py.

◆ dataIds

lsst.ip.diffim.getTemplate.dataIds = collections.defaultdict(list)

Definition at line 208 of file getTemplate.py.

◆ detectorPolygon

lsst.ip.diffim.getTemplate.detectorPolygon = geom.Box2D(bbox)

Definition at line 205 of file getTemplate.py.

◆ overlappingArea

int lsst.ip.diffim.getTemplate.overlappingArea = 0

Definition at line 206 of file getTemplate.py.

◆ patchBBox

lsst.ip.diffim.getTemplate.patchBBox = skymap[dataId['tract']][dataId['patch']].getOuterBBox()

Definition at line 213 of file getTemplate.py.

◆ patchCorners

lsst.ip.diffim.getTemplate.patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())

Definition at line 214 of file getTemplate.py.

◆ patchList

lsst.ip.diffim.getTemplate.patchList = dict()

Definition at line 653 of file getTemplate.py.

◆ patchPolygon

lsst.ip.diffim.getTemplate.patchPolygon = afwGeom.Polygon(wcs.skyToPixel(patchCorners))

Definition at line 215 of file getTemplate.py.

◆ patchWcs

lsst.ip.diffim.getTemplate.patchWcs = skymap[dataId['tract']].getWcs()

Definition at line 212 of file getTemplate.py.