Merge pull request #883 from lsst/tickets/DM-42080

DM-42080: Add effTime to computeExposureSummaryStats

python/lsst/pipe/tasks/computeExposureSummaryStats.py

@@ -92,6 +92,32 @@ class ComputeExposureSummaryStatsConfig(pexConfig.Config):
         "robutsness metric calculations (namely, maxDistToNearestPsf and psfTraceRadiusDelta).",
         default=("BAD", "CR", "EDGE", "INTRP", "NO_DATA", "SAT", "SUSPECT"),
     )
+    fiducialSkyBackground = pexConfig.DictField(
+        keytype=str,
+        itemtype=float,
+        doc="Fiducial sky background level (ADU/s) assumed when calculating effective exposure time. "
+        "Keyed by band.",
+        default={'u': 1.0, 'g': 1.0, 'r': 1.0, 'i': 1.0, 'z': 1.0, 'y': 1.0},
+    )
+    fiducialPsfSigma = pexConfig.DictField(
+        keytype=str,
+        itemtype=float,
+        doc="Fiducial PSF sigma (pixels) assumed when calculating effective exposure time. "
+        "Keyed by band.",
+        default={'u': 1.0, 'g': 1.0, 'r': 1.0, 'i': 1.0, 'z': 1.0, 'y': 1.0},
+    )
+    fiducialZeroPoint = pexConfig.DictField(
+        keytype=str,
+        itemtype=float,
+        doc="Fiducial zero point assumed when calculating effective exposure time. "
+        "Keyed by band.",
+        default={'u': 25.0, 'g': 25.0, 'r': 25.0, 'i': 25.0, 'z': 25.0, 'y': 25.0},
+    )
+    maxEffectiveTransparency = pexConfig.Field(
+        dtype=float,
+        doc="Maximum value allowed for effective transparency scale factor (often inf or 1.0).",
+        default=float('inf')
+    )
 
     def setDefaults(self):
         super().setDefaults()
@@ -151,6 +177,12 @@ class ComputeExposureSummaryStatsTask(pipeBase.Task):
     (against, e.g., extrapolation instability):
     - maxDistToNearestPsf
     - psfTraceRadiusDelta
+
+    These quantities are computed to assess depth:
+    - effTime
+    - effTimePsfSigmaScale
+    - effTimeSkyBgScale
+    - effTimeZeroPointScale
     """
     ConfigClass = ComputeExposureSummaryStatsConfig
     _DefaultName = "computeExposureSummaryStats"
@@ -195,6 +227,8 @@ def run(self, exposure, sources, background):
 
         self.update_masked_image_stats(summary, exposure.getMaskedImage())
 
+        self.update_effective_time_stats(summary, exposure)
+
         md = exposure.getMetadata()
         if 'SFM_ASTROM_OFFSET_MEAN' in md:
             summary.astromOffsetMean = md['SFM_ASTROM_OFFSET_MEAN']
@@ -471,6 +505,101 @@ def update_masked_image_stats(self, summary, masked_image):
         meanVar, _ = statObj.getResult(afwMath.MEANCLIP)
         summary.meanVar = meanVar
 
+    def update_effective_time_stats(self, summary, exposure):
+        """Compute effective exposure time statistics to estimate depth.
+
+        The effective exposure time is the equivalent shutter open
+        time that would be needed under nominal conditions to give the
+        same signal-to-noise for a point source as what is achieved by
+        the observation of interest. This metric combines measurements
+        of the point-spread function, the sky brightness, and the
+        transparency.
+
+        .. _teff_definitions:
+
+        The effective exposure time and its subcomponents are defined in [1]_
+
+        References
+        ----------
+
+        .. [1] Neilsen, E.H., Bernstein, G., Gruendl, R., and Kent, S. (2016).
+               Limiting Magnitude, \tau, teff, and Image Quality in DES Year 1
+               https://www.osti.gov/biblio/1250877/
+
+
+        Parameters
+        ----------
+        summary : `lsst.afw.image.ExposureSummaryStats`
+            Summary object to update in-place.
+        exposure : `lsst.afw.image.ExposureF`
+            Exposure to grab band and exposure time metadata
+
+        """
+        self.log.info("Updating effective exposure time")
+
+        nan = float("nan")
+        summary.effTime = nan
+        summary.effTimePsfSigmaScale = nan
+        summary.effTimeSkyBgScale = nan
+        summary.effTimeZeroPointScale = nan
+
+        exposureTime = exposure.getInfo().getVisitInfo().getExposureTime()
+        filterLabel = exposure.getFilter()
+        if (filterLabel is None) or (not filterLabel.hasBandLabel):
+            band = None
+        else:
+            band = filterLabel.bandLabel
+
+        if band is None:
+            self.log.warn("No band associated with exposure; effTime not calculated.")
+            return
+
+        # PSF component
+        if np.isnan(summary.psfSigma):
+            self.log.debug("PSF sigma is NaN")
+            f_eff = nan
+        elif band not in self.config.fiducialPsfSigma:
+            self.log.debug(f"Fiducial PSF value not found for {band}")
+            f_eff = nan
+        else:
+            fiducialPsfSigma = self.config.fiducialPsfSigma[band]
+            f_eff = (summary.psfSigma / fiducialPsfSigma)**-2
+
+        # Transparency component (note that exposure time may be removed from zeropoint)
+        if np.isnan(summary.zeroPoint):
+            self.log.debug("Zero point is NaN")
+            c_eff = nan
+        elif band not in self.config.fiducialZeroPoint:
+            self.log.debug(f"Fiducial zero point value not found for {band}")
+            c_eff = nan
+        else:
+            fiducialZeroPoint = self.config.fiducialZeroPoint[band]
+            zeroPointDiff = fiducialZeroPoint - (summary.zeroPoint - 2.5*np.log10(exposureTime))
+            c_eff = min(10**(-2.0*(zeroPointDiff)/2.5), self.config.maxEffectiveTransparency)
+
+        # Sky brightness component (convert to cts/s)
+        if np.isnan(summary.skyBg):
+            self.log.debug("Sky background is NaN")
+            b_eff = nan
+        elif band not in self.config.fiducialSkyBackground:
+            self.log.debug(f"Fiducial sky background value not found for {band}")
+            b_eff = nan
+        else:
+            fiducialSkyBackground = self.config.fiducialSkyBackground[band]
+            b_eff = fiducialSkyBackground/(summary.skyBg/exposureTime)
+
+        # Effective exposure time scale factor
+        t_eff = f_eff * c_eff * b_eff
+
+        # Effective exposure time (seconds)
+        effectiveTime = t_eff * exposureTime
+
+        # Output quantities
+        summary.effTime = float(effectiveTime)
+        summary.effTimePsfSigmaScale = float(f_eff)
+        summary.effTimeSkyBgScale = float(b_eff)
+        summary.effTimeZeroPointScale = float(c_eff)
+
 
 def maximum_nearest_psf_distance(
     image_mask,

python/lsst/pipe/tasks/postprocess.py

@@ -1570,7 +1570,9 @@ def run(self, visitSummaryRefs):
                              'psfStarDeltaE1Scatter', 'psfStarDeltaE2Scatter',
                              'psfStarDeltaSizeMedian', 'psfStarDeltaSizeScatter',
                              'psfStarScaledDeltaSizeScatter',
-                             'psfTraceRadiusDelta', 'maxDistToNearestPsf']
+                             'psfTraceRadiusDelta', 'maxDistToNearestPsf',
+                             'effTime', 'effTimePsfSigmaScale',
+                             'effTimeSkyBgScale', 'effTimeZeroPointScale']
             ccdEntry = summaryTable[selectColumns].to_pandas().set_index('id')
             # 'visit' is the human readable visit number.
             # 'visitId' is the key to the visitId table. They are the same.

tests/test_computeExposureSummaryStats.py

@@ -43,9 +43,14 @@ def testComputeExposureSummary(self):
         np.random.seed(12345)
 
         # Make an exposure with a noise image
+        band = "i"
+        physical_filter = "test-i"
         exposure = afwImage.ExposureF(100, 100)
+        exposure.setFilter(afwImage.FilterLabel(band=band, physical=physical_filter))
+
+        skyMean = 100.0
         skySigma = 10.0
-        exposure.getImage().getArray()[:, :] = np.random.normal(0.0, skySigma, size=(100, 100))
+        exposure.getImage().getArray()[:, :] = np.random.normal(skyMean, skySigma, size=(100, 100))
         exposure.getVariance().getArray()[:, :] = skySigma**2.
 
         # Set the visitInfo
@@ -86,8 +91,13 @@ def testComputeExposureSummary(self):
         background = afwMath.BackgroundList()
         background.append(backobj)
 
-        # Run the task
+        # Configure and run the task
         expSummaryTask = ComputeExposureSummaryStatsTask()
+        # Configure nominal values for effective time calculation
+        expSummaryTask.config.fiducialZeroPoint = {band: float(zp)}
+        expSummaryTask.config.fiducialPsfSigma = {band: float(psfSize)}
+        expSummaryTask.config.fiducialSkyBackground = {band: float(skyMean)}
+        # Run the task
         summary = expSummaryTask.run(exposure, None, background)
 
         # Test the outputs
@@ -114,12 +124,14 @@ def testComputeExposureSummary(self):
 
         # Need to compare background level and noise
         # These are only approximately 0+/-10 because of the small image
-        self.assertFloatsAlmostEqual(summary.skyBg, -0.079, atol=1e-3)
-
+        self.assertFloatsAlmostEqual(summary.skyBg, skyMean, rtol=1e-3)
         self.assertFloatsAlmostEqual(summary.meanVar, skySigma**2.)
 
         self.assertFloatsAlmostEqual(summary.zenithDistance, 30.57112, atol=1e-5)
 
+        # Effective exposure time
+        self.assertFloatsAlmostEqual(summary.effTime, 1.0, rtol=1e-3)
+
 
 class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
     pass