Skip to content

Commit

Permalink
Move metric from detectAndMeasure to subtractImages
Browse files Browse the repository at this point in the history 
Update tests accordingly
  • Loading branch information
@mrawls
mrawls committed Jun 5, 2024
1 parent 49cefbd commit 1de5f3a
Show file tree
Hide file tree
Showing 4 changed files with 130 additions and 134 deletions.
56 changes: 2 additions & 54 deletions python/lsst/ip/diffim/detectAndMeasure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>.

import numpy as np
from astropy import units as u

import lsst.afw.detection as afwDetection
import lsst.afw.table as afwTable
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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.
Expand Down
56 changes: 56 additions & 0 deletions python/lsst/ip/diffim/subtractImages.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down
81 changes: 1 addition & 80 deletions tests/test_detectAndMeasure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
71 changes: 71 additions & 0 deletions tests/test_subtractTask.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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):

Expand Down

0 comments on commit 1de5f3a

Please sign in to comment.