Add protection against nans and negatives in reference catalog.

python/lsst/meas/astrom/astrometry.py

@@ -24,6 +24,7 @@
 
 import numpy as np
 from astropy import units
+import scipy.stats
 
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
@@ -68,7 +69,8 @@ class AstrometryConfig(RefMatchConfig):
     magnitudeOutlierRejectionNSigma = pexConfig.Field(
         dtype=float,
         doc=("Number of sigma (measured from the distribution) in magnitude "
-             "for a potential match to be rejected during iteration."),
+             "for a potential reference/source match to be rejected during "
+             "iteration."),
         default=3.0,
     )
 
@@ -397,9 +399,9 @@ def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn):
         Parameters
         ----------
         sourceFluxField : `str`
-            Field in source catalog for fluxes.
+            Field in source catalog for instrumental fluxes.
         refFluxField : `str`
-            Field in reference catalog for fluxes.
+            Field in reference catalog for fluxes (nJy).
         matchesIn : `list` [`lsst.afw.table.ReferenceMatch`]
             List of source/reference matches input
 
@@ -417,20 +419,28 @@ def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn):
             refMag[i] = (match[0][refFluxField]*units.nJy).to_value(units.ABmag)
 
         deltaMag = refMag - sourceMag
-        zp = np.median(deltaMag)
+        # Protect against negative fluxes and nans in the reference catalog.
+        goodDelta, = np.where(np.isfinite(deltaMag))
+        zp = np.median(deltaMag[goodDelta])
         # Use median absolute deviation (MAD) for zpSigma
         # Also require a minimum scatter to prevent floating-point errors from
         # rejecting objects in zero-noise tests.
-        zpSigma = np.clip(1.4826*np.median(np.abs(deltaMag - zp)), 1e-3, None)
+        zpSigma = np.clip(scipy.stats.median_abs_deviation(deltaMag[goodDelta], scale='normal'),
+                          1e-3,
+                          None)
 
-        self.log.debug(f"Rough zeropoint from astrometry matches is {zp:.4f} +/- {zpSigma:.4f}.")
+        self.log.info("Rough zeropoint from astrometry matches is %.4f +/- %.4f.",
+                      zp, zpSigma)
 
-        good, = np.where(np.abs(deltaMag - zp) <= self.config.magnitudeOutlierRejectionNSigma*zpSigma)
-        nOutlier = nMatch - good.size
-        self.log.info(f"Removing {nOutlier} of {nMatch} magnitude outliers.")
+        goodStars = goodDelta[(np.abs(deltaMag[goodDelta] - zp)
+                               <= self.config.magnitudeOutlierRejectionNSigma*zpSigma)]
+
+        nOutlier = nMatch - goodStars.size
+        self.log.info("Removed %d magnitude outliers out of %d total astrometry matches.",
+                      nOutlier, nMatch)
 
         matchesOut = []
-        for matchInd in good:
+        for matchInd in goodStars:
             matchesOut.append(matchesIn[matchInd])
 
         return matchesOut

tests/test_astrometryTask.py

@@ -24,6 +24,8 @@
 import math
 import unittest
 
+from astropy import units
+import scipy.stats
 import numpy as np
 
 import lsst.utils.tests
@@ -222,9 +224,59 @@ def makeSourceCat(self, distortedWcs):
             src.set(sourceCentroidKey, refObj.get(refCentroidKey))
             src.set(sourceInstFluxKey, refObj.get(refFluxRKey))
             src.set(sourceInstFluxErrKey, refObj.get(refFluxRKey)/100)
+
         return sourceCat
 
 
+class TestMagnitudeOutliers(lsst.utils.tests.TestCase):
+    def testMagnitudeOutlierRejection(self):
+        """Test rejection of magnitude outliers."""
+        config = AstrometryTask.ConfigClass()
+        config.doMagnitudeOutlierRejection = True
+        config.magnitudeOutlierRejectionNSigma = 4.0
+        solver = AstrometryTask(config=config, refObjLoader=None)
+
+        nTest = 100
+
+        refSchema = lsst.afw.table.SimpleTable.makeMinimalSchema()
+        refSchema.addField('refFlux', 'F')
+        refCat = lsst.afw.table.SimpleCatalog(refSchema)
+        refCat.resize(nTest)
+
+        srcSchema = lsst.afw.table.SourceTable.makeMinimalSchema()
+        srcSchema.addField('srcFlux', 'F')
+        srcCat = lsst.afw.table.SourceCatalog(srcSchema)
+        srcCat.resize(nTest)
+
+        np.random.seed(12345)
+        refMag = np.full(nTest, 20.0)
+        srcMag = np.random.normal(size=nTest, loc=0.0, scale=1.0)
+
+        # Determine the sigma of the random sample
+        zp = np.median(refMag[: -4] - srcMag[: -4])
+        sigma = scipy.stats.median_abs_deviation(srcMag[: -4], scale='normal')
+
+        # Deliberately alter some magnitudes to be outliers.
+        srcMag[-3] = 4.1*sigma + (20.0 - zp)
+        srcMag[-4] = -4.1*sigma + (20.0 - zp)
+
+        refCat['refFlux'] = (refMag*units.ABmag).to_value(units.nJy)
+        srcCat['srcFlux'] = 10.0**(srcMag/(-2.5))
+
+        # Deliberately poison some reference fluxes.
+        refCat['refFlux'][-1] = np.inf
+        refCat['refFlux'][-2] = np.nan
+
+        matchesIn = []
+        for ref, src in zip(refCat, srcCat):
+            matchesIn.append(lsst.afw.table.ReferenceMatch(first=ref, second=src, distance=0.0))
+
+        matchesOut = solver._removeMagnitudeOutliers('srcFlux', 'refFlux', matchesIn)
+
+        # We should lose the 4 outliers we created.
+        self.assertEqual(len(matchesOut), len(matchesIn) - 4)
+
+
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass