diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 72411e05..d1f5cb5f 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -20,7 +20,6 @@ # along with this program. If not, see . import numpy as np -from astropy import units as u import lsst.afw.detection as afwDetection import lsst.afw.table as afwTable @@ -425,7 +424,7 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor if self.config.doForcedMeasurement: self.measureForcedSources(diaSources, science, difference.getWcs()) - self.calculateMetrics(difference, matchedTemplate, science) + self.calculateMetrics(difference) measurementResults = pipeBase.Struct( subtractedMeasuredExposure=difference, @@ -592,7 +591,7 @@ def measureForcedSources(self, diaSources, science, wcs): for diaSource, forcedSource in zip(diaSources, forcedSources): diaSource.assign(forcedSource, mapper) - def calculateMetrics(self, difference, matchedTemplate, science): + def calculateMetrics(self, difference): """Add difference image QA metrics to the Task metadata. This may be used to produce corresponding metrics (see @@ -602,12 +601,6 @@ def calculateMetrics(self, difference, matchedTemplate, science): ---------- difference : `lsst.afw.image.Exposure` The target difference image to calculate metrics for. - matchedTemplate : `lsst.afw.image.Exposure` - The warped and PSF-matched template used to create the difference - image. - science : `lsst.afw.image.Exposure` - Science exposure that the template was subtracted from. - """ mask = difference.mask badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0 @@ -630,51 +623,6 @@ def calculateMetrics(self, difference, matchedTemplate, science): self.metadata.add("%s_mask_fraction"%maskPlane.lower(), -1) self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane) - maglim_science = self._calculateMagLim(science) - fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy) - maglim_template = self._calculateMagLim(matchedTemplate) - if np.isnan(maglim_template): - self.log.info("Cannot evaluate template limiting mag; adopting science limiting mag for diffim") - maglim_diffim = maglim_science - else: - fluxlim_template = (maglim_template*u.ABmag).to_value(u.nJy) - maglim_diffim = (np.sqrt(fluxlim_science**2 + fluxlim_template**2)*u.nJy).to(u.ABmag).value - self.metadata.add("scienceLimitingMagnitude", maglim_science) - self.metadata.add("templateLimitingMagnitude", maglim_template) - self.metadata.add("diffimLimitingMagnitude", maglim_diffim) - - def _calculateMagLim(self, exposure, nsigma=5.0): - """Calculate an exposure's limiting magnitude. - - This method uses the photometric zeropoint together with the - PSF size from the center of the exposure. - - Parameters - ---------- - exposure : `lsst.afw.image.Exposure` - The target exposure to calculate the limiting magnitude for. - nsigma : `float`, optional - The detection threshold in sigma. - - Returns - ------- - maglim : `astropy.units.Quantity` - The limiting magnitude of the exposure, or np.nan. - """ - try: - psf = exposure.getPsf() - psf_shape = psf.computeShape(psf.getAveragePosition()) - except (InvalidParameterError, afwDetection.InvalidPsfError): - self.log.info("Unable to evaluate PSF, setting maglim to nan") - maglim = np.nan - else: - # Get a more accurate area than `psf_shape.getArea()` via moments - psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy()) - zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1) - maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area)) - finally: - return maglim - def _runStreakMasking(self, maskedImage): """Do streak masking at put results into catalog. diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index f68aa024..a198d127 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -21,8 +21,10 @@ import warnings +from astropy import units as u import numpy as np +import lsst.afw.detection as afwDetection import lsst.afw.image import lsst.afw.math import lsst.geom @@ -407,6 +409,20 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None, self.metadata.add("sciencePsfSize", sciencePsfSize) self.metadata.add("templatePsfSize", templatePsfSize) + # Calculate estimated image depths, i.e., limiting magnitudes + maglim_science = self._calculateMagLim(science, fallbackPsfSize=sciencePsfSize) + fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy) + maglim_template = self._calculateMagLim(template, fallbackPsfSize=templatePsfSize) + if np.isnan(maglim_template): + self.log.info("Cannot evaluate template limiting mag; adopting science limiting mag for diffim") + maglim_diffim = maglim_science + else: + fluxlim_template = (maglim_template*u.ABmag).to_value(u.nJy) + maglim_diffim = (np.sqrt(fluxlim_science**2 + fluxlim_template**2)*u.nJy).to(u.ABmag).value + self.metadata.add("scienceLimitingMagnitude", maglim_science) + self.metadata.add("templateLimitingMagnitude", maglim_template) + self.metadata.add("diffimLimitingMagnitude", maglim_diffim) + if self.config.mode == "auto": convolveTemplate = _shapeTest(template, science, @@ -632,6 +648,46 @@ def finalize(self, template, science, difference, kernel, correctedExposure = difference return correctedExposure + def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None): + """Calculate an exposure's limiting magnitude. + + This method uses the photometric zeropoint together with the + PSF size from the average position of the exposure. + + Parameters + ---------- + exposure : `lsst.afw.image.Exposure` + The target exposure to calculate the limiting magnitude for. + nsigma : `float`, optional + The detection threshold in sigma. + fallbackPsfSize : `float`, optional + PSF FWHM to use in the event the exposure PSF cannot be retrieved. + + Returns + ------- + maglim : `astropy.units.Quantity` + The limiting magnitude of the exposure, or np.nan. + """ + try: + psf = exposure.getPsf() + psf_shape = psf.computeShape(psf.getAveragePosition()) + except (lsst.pex.exceptions.InvalidParameterError, afwDetection.InvalidPsfError): + if fallbackPsfSize is not None: + self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize) + psf_area = np.pi*(fallbackPsfSize/2)**2 + zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1) + maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area)) + else: + self.log.info("Unable to evaluate PSF, setting maglim to nan") + maglim = np.nan + else: + # Get a more accurate area than `psf_shape.getArea()` via moments + psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy()) + zeropoint = exposure.getPhotoCalib().instFluxToMagnitude(1) + maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area)) + finally: + return maglim + @staticmethod def _validateExposures(template, science): """Check that the WCS of the two Exposures match, and the template bbox diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 7e7bc0f8..9ef9bfbf 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -21,13 +21,10 @@ import numpy as np import unittest -from astropy import units as u -import lsst.afw.table as afwTable import lsst.geom from lsst.ip.diffim import detectAndMeasure, subtractImages -from lsst.ip.diffim.utils import makeTestImage, CustomCoaddPsf -import lsst.meas.algorithms as measAlg +from lsst.ip.diffim.utils import makeTestImage from lsst.pipe.base import InvalidQuantumError import lsst.utils.tests @@ -563,82 +560,6 @@ def test_mask_streaks(self): streakMaskSet = (outMask & streakMask) > 0 self.assertTrue(np.all(streakMaskSet[20:23, 40:200])) - def test_metadata_metrics(self): - """Verify fields are added to metadata when detection is run, and - that the difference image limiting magnitude is calculated correctly, - both with a "good" and "bad" seeing template. - """ - science, _ = makeTestImage(psfSize=1.8, doApplyCalibration=True) - matchedTemplate_good, _ = makeTestImage(psfSize=2.4, doApplyCalibration=True) - matchedTemplate_bad, _ = makeTestImage(psfSize=9.5, doApplyCalibration=True) - - # We aren't testing subtraction, so just make a very simple diffim - difference_good = science.clone() - difference_good.maskedImage -= matchedTemplate_good.maskedImage - difference_bad = science.clone() - difference_bad.maskedImage -= matchedTemplate_bad.maskedImage - - # The metadata fields are attached to the detectionTask, so we do - # need to run that; run it for both "good" and "bad" seeing templates - detectionTask_good = self._setup_detection() - detectionTask_bad = self._setup_detection() - _ = detectionTask_good.run(science, matchedTemplate_good, difference_good) - _ = detectionTask_bad.run(science, matchedTemplate_bad, difference_bad) - - # Test that the diffim limiting magnitudes are computed correctly - maglim_science = detectionTask_good._calculateMagLim(science) - fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy) - maglim_template_good = detectionTask_good._calculateMagLim(matchedTemplate_good) - fluxlim_template_good = (maglim_template_good*u.ABmag).to_value(u.nJy) - maglim_template_bad = detectionTask_bad._calculateMagLim(matchedTemplate_bad) - fluxlim_template_bad = (maglim_template_bad*u.ABmag).to_value(u.nJy) - - maglim_good = (np.sqrt(fluxlim_science**2 + fluxlim_template_good**2)*u.nJy).to(u.ABmag).value - maglim_bad = (np.sqrt(fluxlim_science**2 + fluxlim_template_bad**2)*u.nJy).to(u.ABmag).value - - self.assertFloatsAlmostEqual(detectionTask_good.metadata['diffimLimitingMagnitude'], - maglim_good, atol=1e-6) - self.assertFloatsAlmostEqual(detectionTask_bad.metadata['diffimLimitingMagnitude'], - maglim_bad, atol=1e-6) - - # Create a template with a PSF that is not defined at the image center. - # First, make an exposure catalog so we can force the template to have - # a bad (off-image) PSF. It must have a record with a weight field - # and a BBox in order to let us set the PSF manually. - matchedTemplate_offimage, _ = makeTestImage() - schema = afwTable.ExposureTable.makeMinimalSchema() - weightKey = schema.addField("weight", type="D", doc="Coadd weight") - exposureCatalog = afwTable.ExposureCatalog(schema) - record = exposureCatalog.addNew() - record.setD(weightKey, 1.0) - record.setBBox(matchedTemplate_offimage.getBBox()) - kernel = measAlg.DoubleGaussianPsf(7, 7, 2.0).getKernel() - psf = measAlg.KernelPsf(kernel, matchedTemplate_offimage.getBBox().getCenter()) - record.setPsf(psf) - record.setWcs(matchedTemplate_offimage.wcs) - custom_offimage_psf = CustomCoaddPsf(exposureCatalog, matchedTemplate_offimage.wcs) - matchedTemplate_offimage.setPsf(custom_offimage_psf) - - # Test that the bad (off-image) PSF template has a nan maglim, - # and therefore the science maglim is assigned to the diffim. - detectionTask_offimage = self._setup_detection() - difference_offimage = science.clone() - difference_offimage.maskedImage -= matchedTemplate_offimage.maskedImage - _ = detectionTask_offimage.run(science, matchedTemplate_offimage, difference_offimage) - self.assertTrue(np.isnan(detectionTask_offimage.metadata['templateLimitingMagnitude'])) - self.assertEqual(detectionTask_offimage.metadata['diffimLimitingMagnitude'], - detectionTask_offimage.metadata['scienceLimitingMagnitude']) - - # Test that several other expected metadata metrics exist - self.assertIn('nGoodPixels', detectionTask_good.metadata) - self.assertIn('nBadPixels', detectionTask_good.metadata) - self.assertIn('nPixelsDetectedPositive', detectionTask_good.metadata) - self.assertIn('nPixelsDetectedNegative', detectionTask_good.metadata) - self.assertIn('nBadPixelsDetectedPositive', detectionTask_good.metadata) - self.assertIn('nBadPixelsDetectedNegative', detectionTask_good.metadata) - self.assertIn('scienceLimitingMagnitude', detectionTask_good.metadata) - self.assertIn('templateLimitingMagnitude', detectionTask_good.metadata) - class DetectAndMeasureScoreTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase): detectionTask = detectAndMeasure.DetectAndMeasureScoreTask diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index e970440a..59e86970 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -21,6 +21,8 @@ import unittest +from astropy import units as u + import lsst.afw.math as afwMath import lsst.afw.table as afwTable import lsst.geom @@ -878,6 +880,75 @@ def test_fake_mask_plane_propagation(self): self.assertEqual(np.sum(inj_masked.astype(int)-science_fake_masked.astype(int)), 0) self.assertEqual(np.sum(injTmplt_masked.astype(int)-template_fake_masked.astype(int)), 0) + def test_metadata_metrics(self): + """Verify fields are added to metadata when subtraction is run, and + that the difference image limiting magnitude is calculated correctly, + both with a "good" and "bad" seeing template. + """ + science, sources = makeTestImage(psfSize=1.8, doApplyCalibration=True) + template_good, _ = makeTestImage(psfSize=2.4, doApplyCalibration=True) + template_bad, _ = makeTestImage(psfSize=9.5, doApplyCalibration=True) + + # The metadata fields are attached to the subtractTask, so we do + # need to run that; run it for both "good" and "bad" seeing templates + config = subtractImages.AlardLuptonSubtractTask.ConfigClass() + subtractTask_good = subtractImages.AlardLuptonSubtractTask(config=config) + _ = subtractTask_good.run(template_good.clone(), science.clone(), sources) + subtractTask_bad = subtractImages.AlardLuptonSubtractTask(config=config) + _ = subtractTask_bad.run(template_bad.clone(), science.clone(), sources) + + # Test that the diffim limiting magnitudes are computed correctly + maglim_science = subtractTask_good._calculateMagLim(science) + fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy) + maglim_template_good = subtractTask_good._calculateMagLim(template_good) + fluxlim_template_good = (maglim_template_good*u.ABmag).to_value(u.nJy) + maglim_template_bad = subtractTask_bad._calculateMagLim(template_bad) + fluxlim_template_bad = (maglim_template_bad*u.ABmag).to_value(u.nJy) + + maglim_good = (np.sqrt(fluxlim_science**2 + fluxlim_template_good**2)*u.nJy).to(u.ABmag).value + maglim_bad = (np.sqrt(fluxlim_science**2 + fluxlim_template_bad**2)*u.nJy).to(u.ABmag).value + + self.assertFloatsAlmostEqual(subtractTask_good.metadata['diffimLimitingMagnitude'], + maglim_good, atol=1e-6) + self.assertFloatsAlmostEqual(subtractTask_bad.metadata['diffimLimitingMagnitude'], + maglim_bad, atol=1e-6) + + # Create a template with a PSF that is not defined at the image center. + # First, make an exposure catalog so we can force the template to have + # a bad (off-image) PSF. It must have a record with a weight field + # and a BBox in order to let us set the PSF manually. + template_offimage, _ = makeTestImage() + schema = afwTable.ExposureTable.makeMinimalSchema() + weightKey = schema.addField("weight", type="D", doc="Coadd weight") + exposureCatalog = afwTable.ExposureCatalog(schema) + record = exposureCatalog.addNew() + record.setD(weightKey, 1.0) + record.setBBox(template_offimage.getBBox()) + kernel = measAlg.DoubleGaussianPsf(7, 7, 2.0).getKernel() + psf = measAlg.KernelPsf(kernel, template_offimage.getBBox().getCenter()) + record.setPsf(psf) + record.setWcs(template_offimage.wcs) + custom_offimage_psf = CustomCoaddPsf(exposureCatalog, template_offimage.wcs) + template_offimage.setPsf(custom_offimage_psf) + + # Test that template PSF size edge cases are handled correctly. + subtractTask_offimage = subtractImages.AlardLuptonSubtractTask(config=config) + _ = subtractTask_offimage.run(template_offimage.clone(), science.clone(), sources) + # Test that providing no fallbackPsfSize results in a nan template + # limiting magnitude. + maglim_template_offimage = subtractTask_offimage._calculateMagLim(template_offimage) + self.assertTrue(np.isnan(maglim_template_offimage)) + # Test that given the provided fallbackPsfSize, the diffim limiting + # magnitude is calculated correctly. + maglim_template_offimage = 28.3173123103493 + fluxlim_template_offimage = (maglim_template_offimage*u.ABmag).to_value(u.nJy) + maglim_offimage = (np.sqrt(fluxlim_science**2 + fluxlim_template_offimage**2)*u.nJy).to(u.ABmag).value + self.assertEqual(subtractTask_offimage.metadata['diffimLimitingMagnitude'], maglim_offimage) + + # Test that several other expected metadata metrics exist + self.assertIn('scienceLimitingMagnitude', subtractTask_good.metadata) + self.assertIn('templateLimitingMagnitude', subtractTask_good.metadata) + class AlardLuptonPreconvolveSubtractTest(lsst.utils.tests.TestCase):