Merge pull request #449 from lsst/tickets/DM-27896

DM-27896: Add summary statistics component to Exposure

doc/lsst.pipe.tasks/tasks/lsst.pipe.tasks.computeExposureSummary.ComputeExposureSummaryTask.rst

@@ -0,0 +1,26 @@
+.. lsst-task-topic:: lsst.pipe.tasks.computeExposureSummary.ComputeExposureSummaryTask
+
+##########################
+ComputeExposureSummaryTask
+##########################
+
+.. _lsst.pipe.tasks.computeExposureSummary.ComputeExposureSummaryTask-api:
+
+Python API summary
+==================
+
+.. lsst-task-api-summary:: lsst.pipe.tasks.computeExposureSummary.ComputeExposureSummaryTask
+
+.. _lsst.pipe.tasks.computeExposureSummary.ComputeExposureSummaryTask-subtasks:
+
+Retargetable subtasks
+=====================
+
+.. lsst-task-config-subtasks:: lsst.pipe.tasks.computeExposureSummary.ComputeExposureSummaryTask
+
+.. _lsst.pipe.tasks.computeExposureSummary.ComputeExposureSummaryTask-configs:
+
+Configuration fields
+====================
+
+.. lsst-task-config-fields:: lsst.pipe.tasks.computeExposureSummary.ComputeExposureSummaryTask

python/lsst/pipe/tasks/characterizeImage.py

@@ -38,6 +38,7 @@
 from lsst.meas.deblender import SourceDeblendTask
 from .measurePsf import MeasurePsfTask
 from .repair import RepairTask
+from .computeExposureSummaryStats import ComputeExposureSummaryStatsTask
 from lsst.pex.exceptions import LengthError
 
 __all__ = ["CharacterizeImageConfig", "CharacterizeImageTask"]
@@ -157,6 +158,15 @@ class CharacterizeImageConfig(pipeBase.PipelineTaskConfig,
         target=CatalogCalculationTask,
         doc="Subtask to run catalogCalculation plugins on catalog"
     )
+    doComputeSummaryStats = pexConfig.Field(
+        dtype=bool,
+        default=True,
+        doc="Run subtask to measure exposure summary statistics"
+    )
+    computeSummaryStats = pexConfig.ConfigurableField(
+        target=ComputeExposureSummaryStatsTask,
+        doc="Subtask to run computeSummaryStats on exposure"
+    )
     useSimplePsf = pexConfig.Field(
         dtype=bool,
         default=True,
@@ -368,6 +378,8 @@ def __init__(self, butler=None, refObjLoader=None, schema=None, **kwargs):
             self.makeSubtask('measureApCorr', schema=self.schema)
             self.makeSubtask('applyApCorr', schema=self.schema)
         self.makeSubtask('catalogCalculation', schema=self.schema)
+        if self.config.doComputeSummaryStats:
+            self.makeSubtask('computeSummaryStats')
         self._initialFrame = getDebugFrame(self._display, "frame") or 1
         self._frame = self._initialFrame
         self.schema.checkUnits(parse_strict=self.config.checkUnitsParseStrict)
@@ -502,6 +514,11 @@ def run(self, exposure, exposureIdInfo=None, background=None):
             dmeRes.exposure.getInfo().setApCorrMap(apCorrMap)
             self.applyApCorr.run(catalog=dmeRes.sourceCat, apCorrMap=exposure.getInfo().getApCorrMap())
         self.catalogCalculation.run(dmeRes.sourceCat)
+        if self.config.doComputeSummaryStats:
+            summary = self.computeSummaryStats.run(exposure=dmeRes.exposure,
+                                                   sources=dmeRes.sourceCat,
+                                                   background=dmeRes.background)
+            dmeRes.exposure.getInfo().setSummaryStats(summary)
 
         self.display("measure", exposure=dmeRes.exposure, sourceCat=dmeRes.sourceCat)
 

python/lsst/pipe/tasks/computeExposureSummaryStats.py

@@ -0,0 +1,198 @@
+# This file is part of pipe_tasks.
+#
+# 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 warnings
+import numpy as np
+import astropy.units as units
+from astropy.time import Time
+from astropy.coordinates import AltAz, SkyCoord, EarthLocation
+from lsst.daf.base import DateTime
+
+import lsst.pipe.base as pipeBase
+import lsst.pex.config as pexConfig
+import lsst.afw.math as afwMath
+import lsst.afw.image as afwImage
+import lsst.geom
+
+
+__all__ = ("ComputeExposureSummaryStatsTask", "ComputeExposureSummaryStatsConfig")
+
+
+class ComputeExposureSummaryStatsConfig(pexConfig.Config):
+    """Config for ComputeExposureSummaryTask"""
+    sigmaClip = pexConfig.Field(
+        dtype=float,
+        doc="Sigma for outlier rejection for sky noise.",
+        default=3.0,
+    )
+    clipIter = pexConfig.Field(
+        dtype=int,
+        doc="Number of iterations of outlier rejection for sky noise.",
+        default=2,
+    )
+    badMaskPlanes = pexConfig.ListField(
+        dtype=str,
+        doc="Mask planes that, if set, the associated pixel should not be included sky noise calculation.",
+        default=("NO_DATA", "SUSPECT"),
+    )
+
+
+class ComputeExposureSummaryStatsTask(pipeBase.Task):
+    """Task to compute exposure summary statistics.
+
+    This task computes various quantities suitable for DPDD and other
+    downstream processing at the detector centers, including:
+    - psfSigma
+    - psfArea
+    - psfIxx
+    - psfIyy
+    - psfIxy
+    - ra
+    - decl
+    - zenithDistance
+    - zeroPoint
+    - skyBg
+    - skyNoise
+    - meanVar
+    - raCorners
+    - decCorners
+    """
+    ConfigClass = ComputeExposureSummaryStatsConfig
+    _DefaultName = "computeExposureSummaryStats"
+
+    @pipeBase.timeMethod
+    def run(self, exposure, sources, background):
+        """Measure exposure statistics from the exposure, sources, and background.
+
+        Parameters
+        ----------
+        exposure : `lsst.afw.image.ExposureF`
+        sources : `lsst.afw.table.SourceCatalog`
+        background : `lsst.afw.math.BackgroundList`
+
+        Returns
+        -------
+        summary : `lsst.afw.image.ExposureSummary`
+        """
+        self.log.info("Measuring exposure statistics")
+
+        bbox = exposure.getBBox()
+
+        psf = exposure.getPsf()
+        if psf is not None:
+            shape = psf.computeShape(bbox.getCenter())
+            psfSigma = shape.getDeterminantRadius()
+            psfIxx = shape.getIxx()
+            psfIyy = shape.getIyy()
+            psfIxy = shape.getIxy()
+            im = psf.computeKernelImage(bbox.getCenter())
+            # The calculation of effective psf area is taken from
+            # meas_base/src/PsfFlux.cc#L112. See
+            # https://github.com/lsst/meas_base/blob/
+            # 750bffe6620e565bda731add1509507f5c40c8bb/src/PsfFlux.cc#L112
+            psfArea = np.sum(im.array)/np.sum(im.array**2.)
+        else:
+            psfSigma = np.nan
+            psfIxx = np.nan
+            psfIyy = np.nan
+            psfIxy = np.nan
+            psfArea = np.nan
+
+        wcs = exposure.getWcs()
+        if wcs is not None:
+            sph_pts = wcs.pixelToSky(lsst.geom.Box2D(bbox).getCorners())
+            raCorners = [float(sph.getRa().asDegrees()) for sph in sph_pts]
+            decCorners = [float(sph.getDec().asDegrees()) for sph in sph_pts]
+
+            sph_pt = wcs.pixelToSky(bbox.getCenter())
+            ra = sph_pt.getRa().asDegrees()
+            decl = sph_pt.getDec().asDegrees()
+        else:
+            raCorners = [float(np.nan)]*4
+            decCorners = [float(np.nan)]*4
+            ra = np.nan
+            decl = np.nan
+
+        photoCalib = exposure.getPhotoCalib()
+        if photoCalib is not None:
+            zeroPoint = 2.5*np.log10(photoCalib.getInstFluxAtZeroMagnitude())
+        else:
+            zeroPoint = np.nan
+
+        visitInfo = exposure.getInfo().getVisitInfo()
+        date = visitInfo.getDate()
+
+        if date.isValid():
+            # We compute the zenithDistance at the center of the detector rather
+            # than use the boresight value available via the visitInfo, because
+            # the zenithDistance may vary significantly over a large field of view.
+            observatory = visitInfo.getObservatory()
+            loc = EarthLocation(lat=observatory.getLatitude().asDegrees()*units.deg,
+                                lon=observatory.getLongitude().asDegrees()*units.deg,
+                                height=observatory.getElevation()*units.m)
+            obstime = Time(visitInfo.getDate().get(system=DateTime.MJD),
+                           location=loc, format='mjd')
+            coord = SkyCoord(ra*units.degree, decl*units.degree, obstime=obstime, location=loc)
+            with warnings.catch_warnings():
+                warnings.simplefilter('ignore')
+                altaz = coord.transform_to(AltAz)
+
+            zenithDistance = altaz.alt.degree
+        else:
+            zenithDistance = np.nan
+
+        if background is not None:
+            bgStats = (bg[0].getStatsImage().getImage().array
+                       for bg in background)
+            skyBg = sum(np.median(bg[np.isfinite(bg)]) for bg in bgStats)
+        else:
+            skyBg = np.nan
+
+        statsCtrl = afwMath.StatisticsControl()
+        statsCtrl.setNumSigmaClip(self.config.sigmaClip)
+        statsCtrl.setNumIter(self.config.clipIter)
+        statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes))
+        statsCtrl.setNanSafe(True)
+
+        statObj = afwMath.makeStatistics(exposure.getMaskedImage(), afwMath.STDEVCLIP,
+                                         statsCtrl)
+        skyNoise, _ = statObj.getResult(afwMath.STDEVCLIP)
+
+        statObj = afwMath.makeStatistics(exposure.getMaskedImage().getVariance(),
+                                         exposure.getMaskedImage().getMask(),
+                                         afwMath.MEANCLIP, statsCtrl)
+        meanVar, _ = statObj.getResult(afwMath.MEANCLIP)
+
+        summary = afwImage.ExposureSummaryStats(psfSigma=float(psfSigma),
+                                                psfArea=float(psfArea),
+                                                psfIxx=float(psfIxx),
+                                                psfIyy=float(psfIyy),
+                                                psfIxy=float(psfIxy),
+                                                ra=float(ra),
+                                                decl=float(decl),
+                                                zenithDistance=float(zenithDistance),
+                                                zeroPoint=float(zeroPoint),
+                                                skyBg=float(skyBg),
+                                                skyNoise=float(skyNoise),
+                                                meanVar=float(meanVar),
+                                                raCorners=raCorners,
+                                                decCorners=decCorners)
+
+        return summary