Merge pull request #328 from lsst/tickets/DM-43849

DM-43849: Create spatiallySampledMetric to visualize the diffim kernel

python/lsst/ip/diffim/computeSpatiallySampledMetrics.py

@@ -27,7 +27,7 @@
 import lsst.pipe.base as pipeBase
 import lsst.pex.config as pexConfig
 
-from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction
+from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction, getKernelCenterDisplacement
 from lsst.meas.algorithms import SkyObjectsTask
 from lsst.pex.exceptions import InvalidParameterError, RangeError
 from lsst.utils.timer import timeMethod
@@ -72,6 +72,12 @@ class SpatiallySampledMetricsConnections(pipeBase.PipelineTaskConnections,
         storageClass="SourceCatalog",
         name="{fakesType}{coaddName}Diff_candidateDiaSrc",
     )
+    psfMatchingKernel = pipeBase.connectionTypes.Input(
+        doc="Kernel used to PSF match the science and template images.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="MatchingKernel",
+        name="{fakesType}{coaddName}Diff_psfMatchKernel",
+    )
     spatiallySampledMetrics = pipeBase.connectionTypes.Output(
         doc="Summary metrics computed at randomized locations.",
         dimensions=("instrument", "visit", "detector"),
@@ -164,9 +170,32 @@ def __init__(self, **kwargs):
                 "%s_mask_fraction"%maskPlane.lower(), "F",
                 "Fraction of pixels with %s mask"%maskPlane
             )
+        self.schema.addField(
+            "psfMatchingKernel_sum", "F",
+            "PSF matching kernel sum at location.")
+        self.schema.addField(
+            "psfMatchingKernel_dx", "F",
+            "PSF matching kernel centroid offset in x at location.",
+            units="pixel")
+        self.schema.addField(
+            "psfMatchingKernel_dy", "F",
+            "PSF matching kernel centroid offset in y at location.",
+            units="pixel")
+        self.schema.addField(
+            "psfMatchingKernel_length", "F",
+            "PSF matching kernel centroid offset module.",
+            units="arcsecond")
+        self.schema.addField(
+            "psfMatchingKernel_position_angle", "F",
+            "PSF matching kernel centroid offset position angle.",
+            units="radian")
+        self.schema.addField(
+            "psfMatchingKernel_direction", "F",
+            "PSF matching kernel centroid offset direction in detector plane.",
+            units="radian")
 
     @timeMethod
-    def run(self, science, matchedTemplate, template, difference, diaSources):
+    def run(self, science, matchedTemplate, template, difference, diaSources, psfMatchingKernel):
         """Calculate difference image metrics on specific locations across the images
 
         Parameters
@@ -183,6 +212,8 @@ def run(self, science, matchedTemplate, template, difference, diaSources):
             Result of subtracting template from the science image.
         diaSources : `lsst.afw.table.SourceCatalog`
                 The catalog of detected sources.
+        psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`
+            The PSF matching kernel of the subtraction to evaluate.
 
         Returns
         -------
@@ -207,12 +238,13 @@ def run(self, science, matchedTemplate, template, difference, diaSources):
 
         for src in spatiallySampledMetrics:
             self._evaluateLocalMetric(src, science, matchedTemplate, template, difference, diaSources,
-                                      metricsMaskPlanes=metricsMaskPlanes)
+                                      metricsMaskPlanes=metricsMaskPlanes,
+                                      psfMatchingKernel=psfMatchingKernel)
 
         return pipeBase.Struct(spatiallySampledMetrics=spatiallySampledMetrics.asAstropy())
 
     def _evaluateLocalMetric(self, src, science, matchedTemplate, template, difference, diaSources,
-                             metricsMaskPlanes):
+                             metricsMaskPlanes, psfMatchingKernel):
         """Calculate image quality metrics at spatially sampled locations.
 
         Parameters
@@ -229,6 +261,8 @@ def _evaluateLocalMetric(self, src, science, matchedTemplate, template, differen
             Result of subtracting template from the science image.
         metricsMaskPlanes : `list` of `str`
             Mask planes to calculate metrics from.
+        psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`
+            The PSF matching kernel of the subtraction to evaluate.
         """
         bbox = src.getFootprint().getBBox()
         pix = bbox.getCenter()
@@ -271,3 +305,24 @@ def _evaluateLocalMetric(self, src, science, matchedTemplate, template, differen
             src.set("%s_mask_fraction"%maskPlane.lower(),
                     evaluateMaskFraction(difference.mask[bbox], maskPlane)
                     )
+
+        krnlSum, dx, dy, direction, length = getKernelCenterDisplacement(
+            psfMatchingKernel, src.get('x'), src.get('y'))
+
+        point1 = lsst.geom.SpherePoint(
+            src.get('coord_ra'), src.get('coord_dec'),
+            lsst.geom.radians)
+        point2 = science.wcs.pixelToSky(src.get('x') + dx, src.get('y') + dy)
+        bearing = point1.bearingTo(point2)
+        pa_ref_angle = lsst.geom.Angle(np.pi/2, lsst.geom.radians)
+        pa = pa_ref_angle - bearing
+        # Wrap around to get Delta_RA from -pi to +pi
+        pa = pa.wrapCtr()
+        position_angle = pa.asRadians()
+
+        src.set('psfMatchingKernel_sum', krnlSum)
+        src.set('psfMatchingKernel_dx', dx)
+        src.set('psfMatchingKernel_dy', dy)
+        src.set('psfMatchingKernel_length', length*pixScale.asArcseconds())
+        src.set('psfMatchingKernel_position_angle', position_angle)  # in E of N position angle
+        src.set('psfMatchingKernel_direction', direction)  # direction offset in detector

python/lsst/ip/diffim/utils.py

@@ -22,7 +22,7 @@
 """Support utilities for Measuring sources"""
 
 # Export DipoleTestImage to expose fake image generating funcs
-__all__ = ["DipoleTestImage", "evaluateMeanPsfFwhm", "getPsfFwhm"]
+__all__ = ["DipoleTestImage", "evaluateMeanPsfFwhm", "getPsfFwhm", "getKernelCenterDisplacement"]
 
 import itertools
 import numpy as np
@@ -853,6 +853,65 @@ def detectDipoleSources(self, doMerge=True, diffim=None, detectSigma=5.5, grow=3
             return detectTask, schema
 
 
+def getKernelCenterDisplacement(kernel, x, y, image=None):
+    """Calculate the PSF matching kernel peak offset from the nominal
+    position.
+
+    Parameters
+    ----------
+    kernel : `~lsst.afw.math.LinearCombinationKernel`
+        The PSF matching kernel to evaluate.
+    x : `float`
+        The x position on the detector to evaluate the kernel
+    y : `float`
+        The y position on the detector to evaluate the kernel
+    image : `~lsst.afw.image.ImageD`
+        The image to use as base for computing kernel pixel values
+
+    Returns
+    -------
+    kernel_sum : `float`
+        The sum of the kernel on the desired location
+    dx : `float`
+        The displacement of the kernel averaged peak, with respect to the
+        center of the extraction of the kernel
+    dy : `float`
+        The displacement of the kernel averaged peak, with respect to the
+        center of the extraction of the kernel
+    pos_angle: `float`
+        The position angle in detector coordinates of the displacement
+    length : `float`
+        The displacement module of the kernel centroid in pixel units
+    """
+
+    if image is None:
+        image = afwImage.ImageD(kernel.getDimensions())
+
+    # obtain the kernel image
+    hsize = kernel.getWidth()//2
+    kernel_sum = kernel.computeImage(image, doNormalize=False, x=x, y=y)
+
+    data = image.array
+    h, w = data.shape
+    xx = np.arange(w)
+    yy = np.arange(h)
+
+    # create sum vectors and estimate weighted average
+    vx = data.sum(axis=0)
+    vx /= vx.sum()
+    dx = np.dot(vx, xx) - hsize
+
+    vy = data.sum(axis=1)
+    vy /= vy.sum()
+    dy = np.dot(vy, yy) - hsize
+
+    # obtain position angle and norm of displacement
+    pos_angle = np.arctan2(dy, dx)
+    length = np.sqrt(dx**2 + dy**2)
+
+    return kernel_sum, dx, dy, pos_angle, length
+
+
 def getPsfFwhm(psf, average=True, position=None):
     """Directly calculate the horizontal and vertical widths
     of a PSF at half its maximum value.