DM-43848: Rewrite the spatially sampled metrics into a standalone pipelineTask.

python/lsst/ip/diffim/__init__.py

@@ -25,6 +25,7 @@
 from .diffimLib import *
 
 # Python code
+from .computeSpatiallySampledMetrics import *
 from .dcrModel import *
 from .psfMatch import *
 from .modelPsfMatch import *

python/lsst/ip/diffim/computeSpatiallySampledMetrics.py

@@ -0,0 +1,273 @@
+# This file is part of ip_diffim.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import numpy as np
+
+import lsst.geom
+
+import lsst.afw.table as afwTable
+import lsst.pipe.base as pipeBase
+import lsst.pex.config as pexConfig
+
+from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction
+from lsst.meas.algorithms import SkyObjectsTask
+from lsst.pex.exceptions import InvalidParameterError
+from lsst.utils.timer import timeMethod
+
+import lsst.utils
+
+__all__ = ["SpatiallySampledMetricsConfig", "SpatiallySampledMetricsTask"]
+
+
+class SpatiallySampledMetricsConnections(pipeBase.PipelineTaskConnections,
+                                         dimensions=("instrument", "visit", "detector"),
+                                         defaultTemplates={"coaddName": "deep",
+                                                           "warpTypeSuffix": "",
+                                                           "fakesType": ""}):
+    science = pipeBase.connectionTypes.Input(
+        doc="Input science exposure.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="ExposureF",
+        name="{fakesType}calexp"
+    )
+    matchedTemplate = pipeBase.connectionTypes.Input(
+        doc="Warped and PSF-matched template used to create the difference image.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="ExposureF",
+        name="{fakesType}{coaddName}Diff_matchedExp",
+    )
+    template = pipeBase.connectionTypes.Input(
+        doc="Warped and not PSF-matched template used to create the difference image.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="ExposureF",
+        name="{fakesType}{coaddName}Diff_templateExp",
+    )
+    difference = pipeBase.connectionTypes.Input(
+        doc="Difference image with detection mask plane filled in.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="ExposureF",
+        name="{fakesType}{coaddName}Diff_differenceExp",
+    )
+    diaSources = pipeBase.connectionTypes.Input(
+        doc="Filtered diaSources on the difference image.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="SourceCatalog",
+        name="{fakesType}{coaddName}Diff_candidateDiaSrc",
+    )
+    spatiallySampledMetrics = pipeBase.connectionTypes.Output(
+        doc="Summary metrics computed at randomized locations.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="ArrowAstropy",
+        name="{fakesType}{coaddName}Diff_spatiallySampledMetrics",
+    )
+
+
+class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig,
+                                    pipelineConnections=SpatiallySampledMetricsConnections):
+    """Config for SpatiallySampledMetricsTask
+    """
+    metricsMaskPlanes = lsst.pex.config.ListField(
+        dtype=str,
+        doc="List of mask planes to include in metrics",
+        default=('BAD', 'CLIPPED', 'CR', 'DETECTED', 'DETECTED_NEGATIVE', 'EDGE',
+                 'INEXACT_PSF', 'INJECTED', 'INJECTED_TEMPLATE', 'INTRP', 'NOT_DEBLENDED',
+                 'NO_DATA', 'REJECTED', 'SAT', 'SAT_TEMPLATE', 'SENSOR_EDGE', 'STREAK', 'SUSPECT',
+                 'UNMASKEDNAN',
+                 ),
+    )
+    metricSources = pexConfig.ConfigurableField(
+        target=SkyObjectsTask,
+        doc="Generate QA metric sources",
+    )
+
+    def setDefaults(self):
+        self.metricSources.avoidMask = ["NO_DATA", "EDGE"]
+
+
+class SpatiallySampledMetricsTask(lsst.pipe.base.PipelineTask):
+    """Detect and measure sources on a difference image.
+    """
+    ConfigClass = SpatiallySampledMetricsConfig
+    _DefaultName = "spatiallySampledMetrics"
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self.makeSubtask("metricSources")
+        self.schema = afwTable.SourceTable.makeMinimalSchema()
+        self.schema.addField(
+            "x", "F",
+            "X location of the metric evaluation.",
+            units="pixel")
+        self.schema.addField(
+            "y", "F",
+            "Y location of the metric evaluation.",
+            units="pixel")
+        self.metricSources.skySourceKey = self.schema.addField("sky_source", type="Flag",
+                                                               doc="Metric evaluation objects.")
+        self.schema.addField(
+            "source_density", "F",
+            "Density of diaSources at location.",
+            units="count/degree^2")
+        self.schema.addField(
+            "dipole_density", "F",
+            "Density of dipoles at location.",
+            units="count/degree^2")
+        self.schema.addField(
+            "dipole_direction", "F",
+            "Mean dipole orientation.",
+            units="radian")
+        self.schema.addField(
+            "dipole_separation", "F",
+            "Mean dipole separation.",
+            units="pixel")
+        self.schema.addField(
+            "template_value", "F",
+            "Median of template at location.",
+            units="nJy")
+        self.schema.addField(
+            "science_value", "F",
+            "Median of science at location.",
+            units="nJy")
+        self.schema.addField(
+            "diffim_value", "F",
+            "Median of diffim at location.",
+            units="nJy")
+        self.schema.addField(
+            "science_psfSize", "F",
+            "Width of the science image PSF at location.",
+            units="pixel")
+        self.schema.addField(
+            "template_psfSize", "F",
+            "Width of the template image PSF at location.",
+            units="pixel")
+        for maskPlane in self.config.metricsMaskPlanes:
+            self.schema.addField(
+                "%s_mask_fraction"%maskPlane.lower(), "F",
+                "Fraction of pixels with %s mask"%maskPlane
+            )
+
+    @timeMethod
+    def run(self, science, matchedTemplate, template, difference, diaSources):
+        """Calculate difference image metrics on specific locations across the images
+
+        Parameters
+        ----------
+        science : `lsst.afw.image.ExposureF`
+            Science exposure that the template was subtracted from.
+        matchedTemplate : `lsst.afw.image.ExposureF`
+            Warped and PSF-matched template that was used produce the
+            difference image.
+        template : `lsst.afw.image.ExposureF`
+            Warped and non PSF-matched template that was used produce
+            the difference image.
+        difference : `lsst.afw.image.ExposureF`
+            Result of subtracting template from the science image.
+        diaSources : `lsst.afw.table.SourceCatalog`
+                The catalog of detected sources.
+
+        Returns
+        -------
+        results : `lsst.pipe.base.Struct`
+            ``spatiallySampledMetrics`` :  `astropy.table.Table`
+                Image quality metrics spatially sampled locations.
+        """
+
+        idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
+
+        spatiallySampledMetrics = afwTable.SourceCatalog(self.schema)
+        spatiallySampledMetrics.getTable().setIdFactory(idFactory)
+
+        self.metricSources.run(mask=science.mask, seed=difference.info.id, catalog=spatiallySampledMetrics)
+
+        metricsMaskPlanes = []
+        for maskPlane in self.config.metricsMaskPlanes:
+            try:
+                metricsMaskPlanes.append(maskPlane)
+            except InvalidParameterError:
+                self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane)
+
+        for src in spatiallySampledMetrics:
+            self._evaluateLocalMetric(src, science, matchedTemplate, template, difference, diaSources,
+                                      metricsMaskPlanes=metricsMaskPlanes)
+
+        return pipeBase.Struct(spatiallySampledMetrics=spatiallySampledMetrics.asAstropy())
+
+    def _evaluateLocalMetric(self, src, science, matchedTemplate, template, difference, diaSources,
+                             metricsMaskPlanes):
+        """Calculate image quality metrics at spatially sampled locations.
+
+        Parameters
+        ----------
+        src : `lsst.afw.table.SourceRecord`
+            The source record to be updated with metric calculations.
+        diaSources : `lsst.afw.table.SourceCatalog`
+            The catalog of detected sources.
+        science : `lsst.afw.image.Exposure`
+            The science image.
+        matchedTemplate : `lsst.afw.image.Exposure`
+            The reference image, warped and psf-matched to the science image.
+        difference : `lsst.afw.image.Exposure`
+            Result of subtracting template from the science image.
+        metricsMaskPlanes : `list` of `str`
+            Mask planes to calculate metrics from.
+        """
+        bbox = src.getFootprint().getBBox()
+        pix = bbox.getCenter()
+        src.set('science_psfSize', getPsfFwhm(science.psf, position=pix))
+        try:
+            src.set('template_psfSize', getPsfFwhm(template.psf, position=pix))
+        except InvalidParameterError:
+            src.set('template_psfSize', np.nan)
+
+        metricRegionSize = 100
+        bbox.grow(metricRegionSize)
+        bbox = bbox.clippedTo(science.getBBox())
+        nPix = bbox.getArea()
+        pixScale = science.wcs.getPixelScale()
+        area = nPix*pixScale.asDegrees()**2
+        peak = src.getFootprint().getPeaks()[0]
+        src.set('x', peak['i_x'])
+        src.set('y', peak['i_y'])
+        src.setCoord(science.wcs.pixelToSky(peak['i_x'], peak['i_y']))
+        selectSources = diaSources[bbox.contains(diaSources.getX(), diaSources.getY())]
+        sourceDensity = len(selectSources)/area
+        dipoleSources = selectSources[selectSources["ip_diffim_DipoleFit_flag_classification"]]
+        dipoleDensity = len(dipoleSources)/area
+
+        if dipoleSources:
+            meanDipoleOrientation = angleMean(dipoleSources["ip_diffim_DipoleFit_orientation"])
+            src.set('dipole_direction', meanDipoleOrientation)
+            meanDipoleSeparation = np.mean(dipoleSources["ip_diffim_DipoleFit_separation"])
+            src.set('dipole_separation', meanDipoleSeparation)
+
+        templateVal = np.median(matchedTemplate[bbox].image.array)
+        scienceVal = np.median(science[bbox].image.array)
+        diffimVal = np.median(difference[bbox].image.array)
+        src.set('source_density', sourceDensity)
+        src.set('dipole_density', dipoleDensity)
+        src.set('template_value', templateVal)
+        src.set('science_value', scienceVal)
+        src.set('diffim_value', diffimVal)
+        for maskPlane in metricsMaskPlanes:
+            src.set("%s_mask_fraction"%maskPlane.lower(),
+                    evaluateMaskFraction(difference.mask[bbox], maskPlane)
+                    )