LSST Applications g042eb84c57+730a74494b,g04e9c324dd+8c5ae1fdc5,g134cb467dc+1f1e3e7524,g199a45376c+0ba108daf9,g1fd858c14a+fa7d31856b,g210f2d0738+f66ac109ec,g262e1987ae+83a3acc0e5,g29ae962dfc+d856a2cb1f,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+a1e0c9f713,g47891489e3+0d594cb711,g4d44eb3520+c57ec8f3ed,g4d7b6aa1c5+f66ac109ec,g53246c7159+8c5ae1fdc5,g56a1a4eaf3+fd7ad03fde,g64539dfbff+f66ac109ec,g67b6fd64d1+0d594cb711,g67fd3c3899+f66ac109ec,g6985122a63+0d594cb711,g74acd417e5+3098891321,g786e29fd12+668abc6043,g81db2e9a8d+98e2ab9f28,g87389fa792+8856018cbb,g89139ef638+0d594cb711,g8d7436a09f+80fda9ce03,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+a8a29bda4b,g99822b682c+e3ec3c61f9,g9d5c6a246b+0d5dac0c3d,ga41d0fce20+9243b26dd2,gbf99507273+8c5ae1fdc5,gd7ef33dd92+0d594cb711,gdab6d2f7ff+3098891321,ge410e46f29+0d594cb711,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.38
LSST Data Management Base Package
Loading...
Searching...
No Matches
lsst.ip.diffim.getTemplate Namespace Reference

Classes

class  GetDcrTemplateConnections
 
class  GetTemplateConnections
 

Functions

 run (self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter)
 
 checkHighVariance (self, template)
 
 _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()
 
 subfilter = dataId["subfilter"]
 

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 452 of file getTemplate.py.

452 def _makeExposureCatalog(self, exposureRefs, dataIds):
453 """Make an exposure catalog for one tract.
454
455 Parameters
456 ----------
457 exposureRefs : `list` of [`lsst.daf.butler.DeferredDatasetHandle` of \
458 `lsst.afw.image.Exposure`]
459 Exposures to include in the catalog.
460 dataIds : `list` [`lsst.daf.butler.DataCoordinate`]
461 Data ids of each of the included exposures; must have "tract" and
462 "patch" entries.
463
464 Returns
465 -------
466 images : `dict` [`lsst.afw.image.MaskedImage`]
467 MaskedImages of each of the input exposures, for warping.
468 catalog : `lsst.afw.table.ExposureCatalog`
469 Catalog of metadata for each exposure
470 totalBox : `lsst.geom.Box2I`
471 The union of the bounding boxes of all the input exposures.
472 """
473 catalog = afwTable.ExposureCatalog(self.schema)
474 catalog.reserve(len(exposureRefs))
475 exposures = (exposureRef.get() for exposureRef in exposureRefs)
476 images = {}
477 totalBox = geom.Box2I()
478
479 for coadd, dataId in zip(exposures, dataIds):
480 images[dataId] = coadd.maskedImage
481 bbox = coadd.getBBox()
482 totalBox = totalBox.expandedTo(bbox)
483 record = catalog.addNew()
484 record.setPsf(coadd.psf)
485 record.setWcs(coadd.wcs)
486 record.setPhotoCalib(coadd.photoCalib)
487 record.setBBox(bbox)
488 record.setValidPolygon(afwGeom.Polygon(geom.Box2D(bbox).getCorners()))
489 record.set("tract", dataId["tract"])
490 record.set("patch", dataId["patch"])
491 # Weight is used by CoaddPsf, but the PSFs from overlapping patches
492 # should be very similar, so this value mostly shouldn't matter.
493 record.set("weight", 1)
494
495 return images, catalog, totalBox
496
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 571 of file getTemplate.py.

571 def _makePsf(self, template, catalog, wcs):
572 """Return a PSF containing the PSF at each of the input regions.
573
574 Note that although this includes all the exposures from the catalog,
575 the PSF knows which part of the template the inputs came from, so when
576 evaluated at a given position it will not include inputs that never
577 went in to those pixels.
578
579 Parameters
580 ----------
581 template : `lsst.afw.image.Exposure`
582 Generated template the PSF is for.
583 catalog : `lsst.afw.table.ExposureCatalog`
584 Catalog of exposures that went into the template that contains all
585 of the input PSFs.
586 wcs : `lsst.afw.geom.SkyWcs`
587 WCS of the template, to warp the PSFs to.
588
589 Returns
590 -------
591 coaddPsf : `lsst.meas.algorithms.CoaddPsf`
592 The meta-psf constructed from all of the input catalogs.
593 """
594 # CoaddPsf centroid not only must overlap image, but must overlap the
595 # part of image with data. Use centroid of region with data.
596 boolmask = template.mask.array & template.mask.getPlaneBitMask("NO_DATA") == 0
597 maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel))
598 centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
599
600 ctrl = self.config.coaddPsf.makeControl()
601 coaddPsf = CoaddPsf(
602 catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize
603 )
604 return coaddPsf
605
606

◆ _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 497 of file getTemplate.py.

497 def _merge(self, maskedImages, bbox, wcs):
498 """Merge the images that came from one tract into one larger image,
499 ignoring NaN pixels and non-finite variance pixels from individual
500 exposures.
501
502 Parameters
503 ----------
504 maskedImages : `dict` [`lsst.afw.image.MaskedImage` or
505 `lsst.afw.image.Exposure`]
506 Images to be merged into one larger bounding box.
507 bbox : `lsst.geom.Box2I`
508 Bounding box defining the image to merge into.
509 wcs : `lsst.afw.geom.SkyWcs`
510 WCS of all of the input images to set on the output image.
511
512 Returns
513 -------
514 merged : `lsst.afw.image.MaskedImage`
515 Merged image with all of the inputs at their respective bbox
516 positions.
517 count : `int`
518 Count of the number of good pixels (those with positive weights)
519 in the merged image.
520 included : `list` [`int`]
521 List of indexes of patches that were included in the merged
522 result, to be used to trim the exposure catalog.
523 """
524 merged = afwImage.ExposureF(bbox, wcs)
525 weights = afwImage.ImageF(bbox)
526 included = [] # which patches were included in the result
527 for i, (dataId, maskedImage) in enumerate(maskedImages.items()):
528 # Only merge into the trimmed box, to save memory
529 clippedBox = geom.Box2I(maskedImage.getBBox())
530 clippedBox.clip(bbox)
531 if clippedBox.area == 0:
532 self.log.debug("%s does not overlap template region.", dataId)
533 continue # nothing in this image overlaps the output
534 maskedImage = maskedImage.subset(clippedBox)
535 # Catch both zero-value and NaN variance plane pixels
536 good = (maskedImage.variance.array > 0) & (
537 np.isfinite(maskedImage.variance.array)
538 )
539 weight = maskedImage.variance.array[good] ** (-0.5)
540 bad = np.isnan(maskedImage.image.array) | ~good
541 # Note that modifying the patch MaskedImage in place is fine;
542 # we're throwing it away at the end anyway.
543 maskedImage.image.array[bad] = 0.0
544 maskedImage.variance.array[bad] = 0.0
545 # Reset mask, too, since these pixels don't contribute to sum.
546 maskedImage.mask.array[bad] = 0
547 # Cannot use `merged.maskedImage *= weight` because that operator
548 # multiplies the variance by the weight twice; in this case
549 # `weight` are the exact values we want to scale by.
550 maskedImage.image.array[good] *= weight
551 maskedImage.variance.array[good] *= weight
552 weights[clippedBox].array[good] += weight
553 # Free memory before creating new large arrays
554 del weight
555 merged.maskedImage[clippedBox] += maskedImage
556 included.append(i)
557
558 good = weights.array > 0
559
560 # Cannot use `merged.maskedImage /= weights` because that
561 # operator divides the variance by the weight twice; in this case
562 # `weights` are the exact values we want to scale by.
563 weights = weights.array[good]
564 merged.image.array[good] /= weights
565 merged.variance.array[good] /= weights
566
567 merged.mask.array[~good] |= merged.mask.getPlaneBitMask("NO_DATA")
568
569 return merged, good.sum(), included
570

◆ _selectDataRef()

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

Definition at line 849 of file getTemplate.py.

849def _selectDataRef(coaddRef, tract, patch):
850 condition = (coaddRef.dataId["tract"] == tract) & (
851 coaddRef.dataId["patch"] == patch
852 )
853 return condition

◆ checkHighVariance()

lsst.ip.diffim.getTemplate.checkHighVariance ( self,
template )
Set a mask plane for regions with unusually high variance.

Parameters
----------
template : `lsst.afw.image.Exposure`
    The warped template exposure, which will be modified in place.

Definition at line 392 of file getTemplate.py.

392 def checkHighVariance(self, template):
393 """Set a mask plane for regions with unusually high variance.
394
395 Parameters
396 ----------
397 template : `lsst.afw.image.Exposure`
398 The warped template exposure, which will be modified in place.
399 """
400 highVarianceMaskPlaneBit = template.mask.addMaskPlane("HIGH_VARIANCE")
401 template.mask.getPlaneBitMask("HIGH_VARIANCE")
402 varianceExposure = template.clone()
403 varianceExposure.image.array = varianceExposure.variance.array
404 varianceBackground = self.varianceBackground.run(varianceExposure).background.getImage().array
405 threshold = self.config.highVarianceThreshold*np.nanmedian(varianceBackground)
406 highVariancePix = varianceBackground > threshold
407 template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit
408

◆ 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 788 of file getTemplate.py.

788 def checkPatchList(self, patchList):
789 """Check that all of the DcrModel subfilters are present for each
790 patch.
791
792 Parameters
793 ----------
794 patchList : `dict`
795 Dict of the patches containing valid data for each tract.
796
797 Raises
798 ------
799 RuntimeError
800 If the number of exposures found for a patch does not match the
801 number of subfilters.
802 """
803 for tract in patchList:
804 for patch in set(patchList[tract]):
805 if patchList[tract].count(patch) != self.config.numSubfilters:
806 raise RuntimeError(
807 "Invalid number of DcrModel subfilters found: %d vs %d expected",
808 patchList[tract].count(patch),
809 self.config.numSubfilters,
810 )
811

◆ 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 812 of file getTemplate.py.

812 def getDcrModel(self, patchList, coaddRefs, visitInfo):
813 """Build DCR-matched coadds from a list of exposure references.
814
815 Parameters
816 ----------
817 patchList : `dict`
818 Dict of the patches containing valid data for each tract.
819 coaddRefs : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
820 Data references to `~lsst.afw.image.Exposure` representing
821 DcrModels that overlap the detector.
822 visitInfo : `lsst.afw.image.VisitInfo`
823 Metadata for the science image.
824
825 Returns
826 -------
827 coaddExposures : `list` [`lsst.afw.image.Exposure`]
828 Coadd exposures that overlap the detector.
829 """
830 coaddExposures = collections.defaultdict(list)
831 for tract in patchList:
832 for patch in set(patchList[tract]):
833 coaddRefList = [
834 coaddRef
835 for coaddRef in coaddRefs
836 if _selectDataRef(coaddRef, tract, patch)
837 ]
838
839 dcrModel = DcrModel.fromQuantum(
840 coaddRefList,
841 self.config.effectiveWavelength,
842 self.config.bandwidth,
843 self.config.numSubfilters,
844 )
845 coaddExposures[tract].append(dcrModel.buildMatchedExposureHandle(visitInfo=visitInfo))
846 return coaddExposures
847
848

◆ 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 280 of file getTemplate.py.

280 def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter):
281 """Warp coadds from multiple tracts and patches to form a template to
282 subtract from a science image.
283
284 Tract and patch overlap regions are combined by a variance-weighted
285 average, and the variance planes are combined with the same weights,
286 not added in quadrature; the overlap regions are not statistically
287 independent, because they're derived from the same original data.
288 The PSF on the template is created by combining the CoaddPsf on each
289 template image into a meta-CoaddPsf.
290
291 Parameters
292 ----------
293 coaddExposureHandles : `dict` [`int`, `list` of \
294 [`lsst.daf.butler.DeferredDatasetHandle` of \
295 `lsst.afw.image.Exposure`]]
296 Coadds to be mosaicked, indexed on tract id.
297 bbox : `lsst.geom.Box2I`
298 Template Bounding box of the detector geometry onto which to
299 resample the ``coaddExposureHandles``. Modified in-place to include the
300 template border.
301 wcs : `lsst.afw.geom.SkyWcs`
302 Template WCS onto which to resample the ``coaddExposureHandles``.
303 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
304 Record of the tract and patch of each coaddExposure, indexed on
305 tract id.
306 physical_filter : `str`
307 Physical filter of the science image.
308
309 Returns
310 -------
311 result : `lsst.pipe.base.Struct`
312 A struct with attributes:
313
314 ``template``
315 A template coadd exposure assembled out of patches
316 (`lsst.afw.image.ExposureF`).
317
318 Raises
319 ------
320 NoWorkFound
321 If no coadds are found with sufficient un-masked pixels.
322 """
323 band, photoCalib = self._checkInputs(dataIds, coaddExposureHandles)
324
325 bbox.grow(self.config.templateBorderSize)
326
327 warped = {}
328 catalogs = []
329 for tract in coaddExposureHandles:
330 maskedImages, catalog, totalBox = self._makeExposureCatalog(
331 coaddExposureHandles[tract], dataIds[tract]
332 )
333 warpedBox = computeWarpedBBox(catalog[0].wcs, bbox, wcs)
334 warpedBox.grow(5) # to ensure we catch all relevant input pixels
335 # Combine images from individual patches together.
336 unwarped, count, included = self._merge(
337 maskedImages, warpedBox, catalog[0].wcs
338 )
339 # Delete `maskedImages` after combining into one large image to reduce peak memory use
340 del maskedImages
341 if count == 0:
342 self.log.info(
343 "No valid pixels from coadd patches in tract %s; not including in output.",
344 tract,
345 )
346 continue
347 warpedBox.clip(totalBox)
348 potentialInput = self.warper.warpExposure(
349 wcs, unwarped.subset(warpedBox), destBBox=bbox
350 )
351
352 # Delete the single large `unwarped` image after warping to reduce peak memory use
353 del unwarped
354 if np.all(
355 potentialInput.mask.array
356 & potentialInput.mask.getPlaneBitMask("NO_DATA")
357 ):
358 self.log.info(
359 "No overlap from coadd patches in tract %s; not including in output.",
360 tract,
361 )
362 continue
363
364 # Trim the exposure catalog to just the patches that were used.
365 tempCatalog = afwTable.ExposureCatalog(self.schema)
366 tempCatalog.reserve(len(included))
367 for i in included:
368 tempCatalog.append(catalog[i])
369 catalogs.append(tempCatalog)
370 warped[tract] = potentialInput.maskedImage
371
372 if len(warped) == 0:
373 raise pipeBase.NoWorkFound("No patches found to overlap science exposure.")
374 # At this point, all entries will be valid, so we can ignore included.
375 template, count, _ = self._merge(warped, bbox, wcs)
376 if count == 0:
377 raise pipeBase.NoWorkFound("No valid pixels in warped template.")
378
379 # Make a single catalog containing all the inputs that were accepted.
380 catalog = afwTable.ExposureCatalog(self.schema)
381 catalog.reserve(sum([len(c) for c in catalogs]))
382 for c in catalogs:
383 catalog.extend(c)
384
385 # Set a mask plane for any regions with exceptionally high variance.
386 self.checkHighVariance(template)
387 template.setPsf(self._makePsf(template, catalog, wcs))
388 template.setFilter(afwImage.FilterLabel(band, physical_filter))
389 template.setPhotoCalib(photoCalib)
390 return pipeBase.Struct(template=template)
391
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 253 of file getTemplate.py.

◆ dataId

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

Definition at line 257 of file getTemplate.py.

◆ dataIds

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

Definition at line 254 of file getTemplate.py.

◆ detectorPolygon

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

Definition at line 251 of file getTemplate.py.

◆ overlappingArea

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

Definition at line 252 of file getTemplate.py.

◆ patchBBox

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

Definition at line 259 of file getTemplate.py.

◆ patchCorners

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

Definition at line 260 of file getTemplate.py.

◆ patchList

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

Definition at line 755 of file getTemplate.py.

◆ patchPolygon

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

Definition at line 261 of file getTemplate.py.

◆ patchWcs

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

Definition at line 258 of file getTemplate.py.

◆ subfilter

lsst.ip.diffim.getTemplate.subfilter = dataId["subfilter"]

Definition at line 758 of file getTemplate.py.