Merge pull request #970 from lsst/tickets/DM-45573

DM-45573: Compute magnitude limit for exposure summary stats

python/lsst/pipe/tasks/computeExposureSummaryStats.py

@@ -36,6 +36,7 @@
 import lsst.geom as geom
 from lsst.meas.algorithms import ScienceSourceSelectorTask
 from lsst.utils.timer import timeMethod
+import lsst.ip.isr as ipIsr
 
 
 class ComputeExposureSummaryStatsConfig(pexConfig.Config):
@@ -131,6 +132,11 @@ class ComputeExposureSummaryStatsConfig(pexConfig.Config):
         doc="Maximum value allowed for effective transparency scale factor (often inf or 1.0).",
         default=float('inf')
     )
+    magLimSnr = pexConfig.Field(
+        dtype=float,
+        doc="Signal-to-noise ratio for computing the magnitude limit depth.",
+        default=5.0
+    )
 
     def setDefaults(self):
         super().setDefaults()
@@ -204,6 +210,7 @@ class ComputeExposureSummaryStatsTask(pipeBase.Task):
     - effTimePsfSigmaScale
     - effTimeSkyBgScale
     - effTimeZeroPointScale
+    - magLim
     """
     ConfigClass = ComputeExposureSummaryStatsConfig
     _DefaultName = "computeExposureSummaryStats"
@@ -254,6 +261,8 @@ def run(self, exposure, sources, background):
 
         self.update_masked_image_stats(summary, exposure.getMaskedImage())
 
+        self.update_magnitude_limit_stats(summary, exposure)
+
         self.update_effective_time_stats(summary, exposure)
 
         md = exposure.getMetadata()
@@ -572,11 +581,12 @@ def update_effective_time_stats(self, summary, exposure):
         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.
+        transparency. It assumes that the observation is
+        sky-background dominated.
 
         .. _teff_definitions:
 
-        The effective exposure time and its subcomponents are defined in [1]_
+        The effective exposure time and its subcomponents are defined in [1]_.
 
         References
         ----------
@@ -585,7 +595,6 @@ def update_effective_time_stats(self, summary, exposure):
                Limiting Magnitude, \tau, teff, and Image Quality in DES Year 1
                https://www.osti.gov/biblio/1250877/
 
-
         Parameters
         ----------
         summary : `lsst.afw.image.ExposureSummaryStats`
@@ -594,8 +603,6 @@ def update_effective_time_stats(self, summary, exposure):
             Exposure to grab band and exposure time metadata
 
         """
-        self.log.info("Updating effective exposure time")
-
         nan = float("nan")
         summary.effTime = nan
         summary.effTimePsfSigmaScale = nan
@@ -624,7 +631,8 @@ def update_effective_time_stats(self, summary, exposure):
             fiducialPsfSigma = self.config.fiducialPsfSigma[band]
             f_eff = (summary.psfSigma / fiducialPsfSigma)**-2
 
-        # Transparency component (note that exposure time may be removed from zeropoint)
+        # Transparency component
+        # Note: Assumes that the zeropoint includes the exposure time
         if np.isnan(summary.zeroPoint):
             self.log.debug("Zero point is NaN")
             c_eff = nan
@@ -659,6 +667,65 @@ def update_effective_time_stats(self, summary, exposure):
         summary.effTimeSkyBgScale = float(b_eff)
         summary.effTimeZeroPointScale = float(c_eff)
 
+    def update_magnitude_limit_stats(self, summary, exposure):
+        """Compute a summary statistic for the point-source magnitude limit at
+        a given signal-to-noise ratio using exposure-level metadata.
+
+        The magnitude limit is calculated at a given SNR from the PSF,
+        sky, zeropoint, and readnoise in accordance with SMTN-002 [1]_,
+        LSE-40 [2]_, and DMTN-296 [3]_.
+
+        References
+        ----------
+
+        .. [1] "Calculating LSST limiting magnitudes and SNR" (SMTN-002)
+        .. [2] "Photon Rates and SNR Calculations" (LSE-40)
+        .. [3] "Calculations of Image and Catalog Depth" (DMTN-296)
+
+
+        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
+        """
+        if exposure.getDetector() is None:
+            summary.magLim = float("nan")
+            return
+
+        # Calculate the average readnoise [e-]
+        readNoiseList = list(ipIsr.getExposureReadNoises(exposure).values())
+        readNoise = np.nanmean(readNoiseList)
+        if np.isnan(readNoise):
+            readNoise = 0.0
+            self.log.warn("Read noise set to NaN! Setting readNoise to 0.0 to proceed.")
+
+        # Calculate the average gain [e-/ADU]
+        gainList = list(ipIsr.getExposureGains(exposure).values())
+        gain = np.nanmean(gainList)
+        if np.isnan(gain):
+            self.log.warn("Gain set to NaN! Setting magLim to NaN.")
+            summary.magLim = float("nan")
+            return
+
+        # Get the image units (default to 'adu' if metadata key absent)
+        md = exposure.getMetadata()
+        if md.get("LSST ISR UNIT", "adu") == "electron":
+            gain = 1.0
+
+        # Convert readNoise to image units
+        readNoise /= gain
+
+        # Calculate the limiting magnitude.
+        # Note 1: Assumes that the image and readnoise have the same units
+        # Note 2: Assumes that the zeropoint includes the exposure time
+        magLim = compute_magnitude_limit(summary.psfArea, summary.skyBg,
+                                         summary.zeroPoint, readNoise,
+                                         gain, self.config.magLimSnr)
+
+        summary.magLim = float(magLim)
+
 
 def maximum_nearest_psf_distance(
     image_mask,
@@ -830,3 +897,65 @@ def compute_ap_corr_sigma_scaled_delta(
         ap_corr_sigma_scaled_delta = (np.max(ap_corr_good) - np.min(ap_corr_good))/psfSigma
 
     return ap_corr_sigma_scaled_delta
+
+
+def compute_magnitude_limit(
+        psfArea,
+        skyBg,
+        zeroPoint,
+        readNoise,
+        gain,
+        snr
+):
+    """Compute the expected point-source magnitude limit at a given
+    signal-to-noise ratio given the exposure-level metadata. Based on
+    the signal-to-noise formula provided in SMTN-002 (see LSE-40 for
+    more details on the calculation).
+
+      SNR = C / sqrt( C/g + (B/g + sigma_inst**2) * neff )
+
+    where C is the counts from the source, B is counts from the (sky)
+    background, sigma_inst is the instrumental (read) noise, neff is
+    the effective size of the PSF, and g is the gain in e-/ADU.  Note
+    that input values of ``skyBg``, ``zeroPoint``, and ``readNoise``
+    should all consistently be in electrons or ADU.
+
+    Parameters
+    ----------
+    psfArea : `float`
+        The effective area of the PSF [pix].
+    skyBg : `float`
+        The sky background counts for the exposure [ADU or e-].
+    zeroPoint : `float`
+        The zeropoint (includes exposure time) [ADU or e-].
+    readNoise : `float`
+        The instrumental read noise for the exposure [ADU or e-].
+    gain : `float`
+        The instrumental gain for the exposure [e-/ADU]. The gain should
+         be 1.0 if the skyBg, zeroPoint, and readNoise are in e-.
+    snr : `float`
+        Signal-to-noise ratio at which magnitude limit is calculated.
+
+    Returns
+    -------
+    magnitude_limit : `float`
+        The expected magnitude limit at the given signal to noise.
+
+    """
+    # Effective number of pixels within the PSF calculated directly
+    # from the PSF model.
+    neff = psfArea
+
+    # Instrumental noise (read noise only)
+    sigma_inst = readNoise
+
+    # Total background counts derived from Eq. 45 of LSE-40
+    background = (skyBg/gain + sigma_inst**2) * neff
+    # Source counts to achieve the desired SNR (from quadratic formula)
+    source = (snr**2)/(2*gain) + np.sqrt((snr**4)/(4*gain) + snr**2 * background)
+
+    # Convert to a magnitude using the zeropoint.
+    # Note: Zeropoint currently includes exposure time
+    magnitude_limit = -2.5*np.log10(source) + zeroPoint
+
+    return magnitude_limit

python/lsst/pipe/tasks/postprocess.py

@@ -1229,7 +1229,8 @@ def run(self, visitSummaryRefs):
                              "psfApFluxDelta", "psfApCorrSigmaScaledDelta",
                              "maxDistToNearestPsf",
                              "effTime", "effTimePsfSigmaScale",
-                             "effTimeSkyBgScale", "effTimeZeroPointScale"]
+                             "effTimeSkyBgScale", "effTimeZeroPointScale",
+                             "magLim"]
             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

@@ -29,10 +29,12 @@
 from lsst.afw.detection import GaussianPsf
 import lsst.afw.image as afwImage
 import lsst.afw.math as afwMath
-from lsst.daf.base import DateTime
+from lsst.daf.base import DateTime, PropertyList
 from lsst.afw.coord import Observatory
 from lsst.afw.geom import makeCdMatrix, makeSkyWcs
+from lsst.afw.cameraGeom.testUtils import DetectorWrapper
 from lsst.pipe.tasks.computeExposureSummaryStats import ComputeExposureSummaryStatsTask
+from lsst.pipe.tasks.computeExposureSummaryStats import compute_magnitude_limit
 
 
 class ComputeExposureSummaryTestCase(lsst.utils.tests.TestCase):
@@ -43,11 +45,22 @@ def testComputeExposureSummary(self):
         np.random.seed(12345)
 
         # Make an exposure with a noise image
+        exposure = afwImage.ExposureF(100, 100)
+
         band = "i"
         physical_filter = "test-i"
-        exposure = afwImage.ExposureF(100, 100)
         exposure.setFilter(afwImage.FilterLabel(band=band, physical=physical_filter))
 
+        readNoise = 5.0
+        detector = DetectorWrapper(numAmps=1).detector
+        metadata = PropertyList()
+        metadata.add("LSST ISR UNIT", "electron")
+        for amp in detector.getAmplifiers():
+            metadata.add(f"LSST ISR READNOISE {amp.getName()}", readNoise)
+            metadata.add(f"LSST ISR GAIN {amp.getName()}", 1.0)
+        exposure.setDetector(detector)
+        exposure.setMetadata(metadata)
+
         skyMean = 100.0
         skySigma = 10.0
         exposure.getImage().getArray()[:, :] = np.random.normal(skyMean, skySigma, size=(100, 100))
@@ -102,10 +115,10 @@ def testComputeExposureSummary(self):
 
         # Configure and run the task
         expSummaryTask = ComputeExposureSummaryStatsTask()
-        # Configure nominal values for effective time calculation
-        expSummaryTask.config.fiducialZeroPoint = {band: float(zp)}
+        # Configure nominal values for effective time calculation (normalized to 1s exposure)
+        expSummaryTask.config.fiducialZeroPoint = {band: float(zp - 2.5*np.log10(expTime))}
         expSummaryTask.config.fiducialPsfSigma = {band: float(psfSize)}
-        expSummaryTask.config.fiducialSkyBackground = {band: float(skyMean)}
+        expSummaryTask.config.fiducialSkyBackground = {band: float(skyMean/expTime)}
         # Run the task
         summary = expSummaryTask.run(exposure, None, background)
 
@@ -145,8 +158,47 @@ def testComputeExposureSummary(self):
 
         self.assertFloatsAlmostEqual(summary.zenithDistance, 30.57112, atol=1e-5)
 
-        # Effective exposure time
-        self.assertFloatsAlmostEqual(summary.effTime, 1.0, rtol=1e-3)
+        # Effective exposure time and depth
+        self.assertFloatsAlmostEqual(summary.effTime, expTime, rtol=1e-3)
+        self.assertFloatsAlmostEqual(summary.magLim, 26.584, rtol=1e-3)
+
+    def testComputeMagnitudeLimit(self):
+        """Test the magnitude limit calculation using fiducials from SMTN-002
+        and syseng_throughputs."""
+
+        # Values from syseng_throughputs notebook assuming 30s exposure
+        # consisting of 2x15s snaps each with readnoise of 9e-
+        fwhm_eff_fid = {'g': 0.87, 'r': 0.83, 'i': 0.80}
+        skycounts_fid = {'g': 463.634122, 'r': 988.626863, 'i': 1588.280513}
+        zeropoint_fid = {'g': 28.508375, 'r': 28.360838, 'i': 28.171396}
+        readnoise_fid = {'g': 12.73, 'r': 12.73, 'i': 12.73}
+        # Assumed values from SMTN-002
+        snr = 5
+        gain = 1.0
+        # Output magnitude limit from syseng_throughputs notebook
+        m5_fid = {'g': 24.90, 'r': 24.48, 'i': 24.10}
+
+        for band in ['g', 'r', 'i']:
+            # Translate into DM quantities
+            psfArea = 2.266 * (fwhm_eff_fid[band] / 0.2)**2
+            skyBg = skycounts_fid[band]
+            zeroPoint = zeropoint_fid[band] + 2.5*np.log10(30)
+            readNoise = readnoise_fid[band]
+
+            # Calculate the M5 values
+            m5 = compute_magnitude_limit(psfArea, skyBg, zeroPoint, readNoise, gain, snr)
+            self.assertFloatsAlmostEqual(m5, m5_fid[band], atol=1e-2)
+
+        # Check that input NaN lead to output NaN
+        nan = float('nan')
+        m5 = compute_magnitude_limit(nan, skyBg, zeroPoint, readNoise, gain, snr)
+        self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True)
+        m5 = compute_magnitude_limit(psfArea, nan, zeroPoint, readNoise, gain, snr)
+        self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True)
+        m5 = compute_magnitude_limit(psfArea, skyBg, nan, readNoise, gain, snr)
+        self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True)
+        m5 = compute_magnitude_limit(psfArea, skyBg, zeroPoint, nan, gain, snr)
+        self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True)
 
 
 class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):