DM-56525: Replace the CoaddPSF of the template_detector with a Gaussian

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,

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,
     )

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
     -------

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.
         """