Merge pull request #144 from lsst/tickets/DM-30490

DM-30490: Add magnitude outlier rejection to AstrometryTask.

python/lsst/meas/astrom/astrometry.py

@@ -22,6 +22,9 @@
 
 __all__ = ["AstrometryConfig", "AstrometryTask"]
 
+import numpy as np
+from astropy import units
+import scipy.stats
 
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
@@ -57,6 +60,19 @@ class AstrometryConfig(RefMatchConfig):
         default=0.001,
         min=0,
     )
+    doMagnitudeOutlierRejection = pexConfig.Field(
+        dtype=bool,
+        doc=("If True then a rough zeropoint will be computed from matched sources "
+             "and outliers will be rejected in the iterations."),
+        default=False,
+    )
+    magnitudeOutlierRejectionNSigma = pexConfig.Field(
+        dtype=float,
+        doc=("Number of sigma (measured from the distribution) in magnitude "
+             "for a potential reference/source match to be rejected during "
+             "iteration."),
+        default=3.0,
+    )
 
     def setDefaults(self):
         # Override the default source selector for astrometry tasks
@@ -342,9 +358,14 @@ def _matchAndFitWcs(self, refCat, sourceCat, goodSourceCat, refFluxField, bbox,
                 title="Initial WCS",
             )
 
+        if self.config.doMagnitudeOutlierRejection:
+            matches = self._removeMagnitudeOutliers(sourceFluxField, refFluxField, matchRes.matches)
+        else:
+            matches = matchRes.matches
+
         self.log.debug("Fitting WCS")
         fitRes = self.wcsFitter.fitWcs(
-            matches=matchRes.matches,
+            matches=matches,
             initWcs=wcs,
             bbox=bbox,
             refCat=refCat,
@@ -358,16 +379,68 @@ def _matchAndFitWcs(self, refCat, sourceCat, goodSourceCat, refFluxField, bbox,
             displayAstrometry(
                 refCat=refCat,
                 sourceCat=matchRes.usableSourceCat,
-                matches=matchRes.matches,
+                matches=matches,
                 exposure=exposure,
                 bbox=bbox,
                 frame=frame + 2,
                 title="Fit TAN-SIP WCS",
             )
 
         return pipeBase.Struct(
-            matches=matchRes.matches,
+            matches=matches,
             wcs=fitWcs,
             scatterOnSky=scatterOnSky,
             match_tolerance=matchRes.match_tolerance,
         )
+
+    def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn):
+        """Remove magnitude outliers, computing a simple zeropoint.
+
+        Parameters
+        ----------
+        sourceFluxField : `str`
+            Field in source catalog for instrumental fluxes.
+        refFluxField : `str`
+            Field in reference catalog for fluxes (nJy).
+        matchesIn : `list` [`lsst.afw.table.ReferenceMatch`]
+            List of source/reference matches input
+
+        Returns
+        -------
+        matchesOut : `list` [`lsst.afw.table.ReferenceMatch`]
+            List of source/reference matches with magnitude
+            outliers removed.
+        """
+        nMatch = len(matchesIn)
+        sourceMag = np.zeros(nMatch)
+        refMag = np.zeros(nMatch)
+        for i, match in enumerate(matchesIn):
+            sourceMag[i] = -2.5*np.log10(match[1][sourceFluxField])
+            refMag[i] = (match[0][refFluxField]*units.nJy).to_value(units.ABmag)
+
+        deltaMag = refMag - sourceMag
+        # 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(scipy.stats.median_abs_deviation(deltaMag[goodDelta], scale='normal'),
+                          1e-3,
+                          None)
+
+        self.log.info("Rough zeropoint from astrometry matches is %.4f +/- %.4f.",
+                      zp, zpSigma)
+
+        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 goodStars:
+            matchesOut.append(matchesIn[matchInd])
+
+        return matchesOut

python/lsst/meas/astrom/pessimistic_pattern_matcher_b_3D.py

@@ -1410,7 +1410,7 @@ def _handshake_match(self, matches_src, matches_ref):
         handshake_mask_array : `numpy.ndarray`, (N,)
            Return the array positions where the two match catalogs agree.
         """
-        handshake_mask_array = np.zeros(len(matches_src), dtype=np.bool)
+        handshake_mask_array = np.zeros(len(matches_src), dtype=bool)
 
         for src_match_idx, match in enumerate(matches_src):
             ref_match_idx = np.searchsorted(matches_ref[:, 1], match[1])

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
@@ -162,6 +164,17 @@ def doTest(self, pixelsToTanPixels, order=3):
         )
         self.assertLess(len(resultsRefSelect.matches), len(results.matches))
 
+        # try again, allowing magnitude outlier rejection.
+        config.doMagnitudeOutlierRejection = True
+        solverMagOutlierRejection = AstrometryTask(config=config, refObjLoader=self.refObjLoader)
+        self.exposure.setWcs(distortedWcs)
+        resultsMagOutlierRejection = solverMagOutlierRejection.run(
+            sourceCat=sourceCat,
+            exposure=self.exposure,
+        )
+        self.assertLess(len(resultsMagOutlierRejection.matches), len(resultsRefSelect.matches))
+        config.doMagnitudeOutlierRejection = False
+
         # try again, but without fitting the WCS, no reference selector
         config.referenceSelector.doUnresolved = False
         config.forceKnownWcs = True
@@ -211,9 +224,67 @@ def makeSourceCat(self, distortedWcs):
             src.set(sourceCentroidKey, refObj.get(refCentroidKey))
             src.set(sourceInstFluxKey, refObj.get(refFluxRKey))
             src.set(sourceInstFluxErrKey, refObj.get(refFluxRKey)/100)
+
+        # Deliberately add some outliers to check that the magnitude
+        # outlier rejection code is being run.
+        sourceCat[sourceInstFluxKey][0: 4] *= 1000.0
+
         return sourceCat
 
 
+class TestMagnitudeOutliers(lsst.utils.tests.TestCase):
+    def testMagnitudeOutlierRejection(self):
+        """Test rejection of magnitude outliers.
+
+        This test only tests the outlier rejection, and not any other
+        part of the matching or astrometry fitter.
+        """
+        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] = (config.magnitudeOutlierRejectionNSigma + 0.1)*sigma + (20.0 - zp)
+        srcMag[-4] = -(config.magnitudeOutlierRejectionNSigma + 0.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