diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 694e8b6b..74e7a758 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -21,8 +21,10 @@ import collections import dataclasses +from astropy.stats import gaussian_sigma_to_fwhm import numpy as np +import lsst.afw.detection as afwDetection import lsst.afw.image as afwImage import lsst.geom as geom import lsst.afw.geom as afwGeom @@ -31,10 +33,12 @@ from lsst.afw.math._warper import computeWarpedBBox import lsst.afw.math as afwMath import lsst.pex.config as pexConfig +import lsst.pex.exceptions import lsst.pipe.base as pipeBase from lsst.skymap import BaseSkyMap from lsst.ip.diffim.dcrModel import DcrModel +from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask, ScaleVarianceTask from lsst.utils.timer import timeMethod @@ -45,6 +49,9 @@ "GetDcrTemplateConfig", ] +_PSF_STAMP_SIGMAS = 5 +"""Width, in N sigma, of the Gaussian approximation to the CoaddPsf.""" + class GetTemplateConnections( pipeBase.PipelineTaskConnections, @@ -122,6 +129,17 @@ class GetTemplateConfig( doc="Configuration for CoaddPsf", dtype=CoaddPsfConfig, ) + fwhmExposureGrid = pexConfig.Field( + doc="Grid size to compute the average PSF FWHM in an exposure (pixels)", + dtype=int, + default=10, + ) + fwhmExposureBuffer = pexConfig.Field( + doc="Fractional buffer margin to be left out of all sides of the image during construction" + " of grid to compute average PSF FWHM in an exposure", + dtype=float, + default=0.05, + ) varianceBackground = pexConfig.ConfigurableField( target=SubtractBackgroundTask, doc="Task to estimate the background variance.", @@ -322,8 +340,8 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi 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. + The PSF of the template is a Gaussian with the same width as the + CoaddPsf evaluated at the center of the detector. Parameters ---------- @@ -354,7 +372,8 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi ``template`` A template coadd exposure assembled out of patches, with - pixels in nJy (`lsst.afw.image.ExposureF`). + pixels in nJy. Its PSF is a Gaussian approximation to the + CoaddPsf of its inputs (`lsst.afw.image.ExposureF`). Raises ------ @@ -456,6 +475,7 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi photoCalib = identity template.setPhotoCalib(photoCalib) template.setPsf(self._makePsf(template, catalog, wcs)) + template.setPsf(self._makeGaussianPsf(template)) # Record the input coadd patches as the template's coadd inputs. coaddInputs = afwImage.CoaddInputs(afwTable.ExposureTable.makeMinimalSchema(), self.schema) @@ -774,6 +794,47 @@ def _makePsf(self, template, catalog, wcs): ) return coaddPsf + def _makeGaussianPsf(self, template): + """Return a Gaussian PSF with the width of the template CoaddPsf. + + Parameters + ---------- + template : `lsst.afw.image.Exposure` + Template with the CoaddPsf to approximate. The width is measured + at the average position of the PSF, or averaged over a grid of + positions if the PSF cannot be evaluated there. + + Returns + ------- + psf : `lsst.afw.detection.GaussianPsf` + Gaussian with the width of the template PSF. + + Raises + ------ + ValueError + Raised if the PSF cannot be evaluated at any grid position. + """ + # The CoaddPsf is not defined where the template has no inputs, which + # can include its average position. + try: + fwhm = getPsfFwhm(template.psf) + except lsst.pex.exceptions.Exception: + # Catch a broad range of exceptions, since some are C++ only + # Catching: + # - lsst::geom::SingularTransformException + # - lsst.pex.exceptions.InvalidParameterError + # - lsst.pex.exceptions.RangeError + self.log.info("Unable to evaluate PSF at the average position. " + "Evaluating PSF on a grid of points.") + fwhm = evaluateMeanPsfFwhm(template, + fwhmExposureBuffer=self.config.fwhmExposureBuffer, + fwhmExposureGrid=self.config.fwhmExposureGrid) + sigma = fwhm/gaussian_sigma_to_fwhm + size = 2*int(np.ceil(_PSF_STAMP_SIGMAS*sigma)) + 1 + self.metadata["templatePsfSigma"] = float(sigma) + self.log.info("Approximated the template PSF with a Gaussian of sigma %.3f pixels.", sigma) + return afwDetection.GaussianPsf(size, size, sigma) + class GetDcrTemplateConnections( GetTemplateConnections, diff --git a/python/lsst/ip/diffim/makeKernel.py b/python/lsst/ip/diffim/makeKernel.py index eeb3459e..3ef370df 100644 --- a/python/lsst/ip/diffim/makeKernel.py +++ b/python/lsst/ip/diffim/makeKernel.py @@ -60,7 +60,7 @@ class MakeKernelConfig(PsfMatchConfig): doc="Initial measurements used to feed stars to kernel fitting", ) fwhmExposureGrid = lsst.pex.config.Field( - doc="Grid size to compute the average PSF FWHM in an exposure", + doc="Grid size to compute the average PSF FWHM in an exposure (pixels)", dtype=int, default=10, ) diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index 5d902cdd..6e3823d8 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -154,7 +154,7 @@ def evaluateMeanPsfFwhm(exposure: afwImage.Exposure, during the construction of the grid to compute mean PSF FWHM in an exposure. fwhmExposureGrid : `int` - Grid size to compute the mean FWHM in an exposure. + Grid size to compute the average PSF FWHM in an exposure (pixels). Returns ------- diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index 7673eb0c..3b407eec 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -22,9 +22,12 @@ import collections import itertools import unittest +import unittest.mock +from astropy.stats import gaussian_sigma_to_fwhm import numpy as np +import lsst.afw.detection import lsst.afw.geom import lsst.afw.image import lsst.afw.math @@ -34,6 +37,7 @@ import lsst.ip.diffim import lsst.meas.algorithms import lsst.meas.base.tests +import lsst.pex.exceptions import lsst.pipe.base as pipeBase import lsst.skymap import lsst.utils.tests @@ -196,10 +200,24 @@ def _checkMetadata(self, template, config, box, wcs, nPsfs): self.assertEqual(template.getXY0(), expectedBox.getMin()) self.assertEqual(template.filter.bandLabel, "a") self.assertEqual(template.filter.physicalLabel, "a_test") - self.assertEqual(template.psf.getComponentCount(), nPsfs) + # The template PSF is a Gaussian with the width of the CoaddPsf of + # its inputs. + self.assertIsInstance(template.psf, lsst.afw.detection.GaussianPsf) + coaddPsf = self._makeCoaddPsf(template, config) + self.assertEqual(coaddPsf.getComponentCount(), nPsfs) + position = coaddPsf.getAveragePosition() + self.assertFloatsAlmostEqual(template.psf.getSigma(), + coaddPsf.computeShape(position).getTraceRadius()) self.assertTrue(template.getInfo().hasCoaddInputs()) self.assertEqual(len(template.getInfo().getCoaddInputs().ccds), nPsfs) + def _makeCoaddPsf(self, template, config): + """Return the CoaddPsf that ``run`` approximated with a Gaussian, + rebuilt from the coadd inputs recorded on the template. + """ + task = lsst.ip.diffim.GetTemplateTask(config=config) + return task._makePsf(template, template.getInfo().getCoaddInputs().ccds, template.wcs) + def _checkPixels(self, template, config, box): """Check that the pixel values in the template are close to the original image, calibrated to nJy. @@ -322,6 +340,35 @@ def _remakePatches(self, photoCalib): for tract_id in range(4): self._makePatches(self.skymap.generateTract(tract_id)) + def testGaussianPsfGridFallback(self): + """If the CoaddPsf cannot be evaluated at its average position, the + Gaussian width is averaged over a grid of positions instead. + """ + box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) + task = lsst.ip.diffim.GetTemplateTask() + result = task.run(coaddExposureHandles={0: self.patches[0]}, + bbox=lsst.geom.Box2I(box), + wcs=self.exposure.wcs, + dataIds={0: self.dataIds[0]}, + physical_filter="a_test") + template = result.template + template.setPsf(self._makeCoaddPsf(template, task.config)) + expected = lsst.ip.diffim.utils.evaluateMeanPsfFwhm( + template, + fwhmExposureBuffer=task.config.fwhmExposureBuffer, + fwhmExposureGrid=task.config.fwhmExposureGrid, + )/gaussian_sigma_to_fwhm + + error = lsst.pex.exceptions.InvalidParameterError("No inputs at the average position.") + with unittest.mock.patch("lsst.ip.diffim.getTemplate.getPsfFwhm", side_effect=error): + with self.assertLogs(task.log.name, level="INFO") as cm: + psf = task._makeGaussianPsf(template) + + self.assertIn("grid of points", "\n".join(cm.output)) + self.assertIsInstance(psf, lsst.afw.detection.GaussianPsf) + self.assertFloatsAlmostEqual(psf.getSigma(), expected) + self.assertFloatsAlmostEqual(task.metadata["templatePsfSigma"], expected) + def testRunNoTemplate(self): """A bounding box that doesn't overlap the patches will raise. """