From bcf184e9c3685dee7fd7cebcd1daabc6b4b16af8 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Thu, 10 Oct 2024 04:47:58 -0700 Subject: [PATCH] Compute magnitude limit for exposure summary stats --- .../pipe/tasks/computeExposureSummaryStats.py | 141 +++++++++++++++++- python/lsst/pipe/tasks/postprocess.py | 3 +- tests/test_computeExposureSummaryStats.py | 66 +++++++- 3 files changed, 196 insertions(+), 14 deletions(-) diff --git a/python/lsst/pipe/tasks/computeExposureSummaryStats.py b/python/lsst/pipe/tasks/computeExposureSummaryStats.py index 3f45271e1..ed9351434 100644 --- a/python/lsst/pipe/tasks/computeExposureSummaryStats.py +++ b/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 diff --git a/python/lsst/pipe/tasks/postprocess.py b/python/lsst/pipe/tasks/postprocess.py index 005bd1675..2418056fa 100644 --- a/python/lsst/pipe/tasks/postprocess.py +++ b/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. diff --git a/tests/test_computeExposureSummaryStats.py b/tests/test_computeExposureSummaryStats.py index 4be9446a1..b2debe300 100644 --- a/tests/test_computeExposureSummaryStats.py +++ b/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):