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,26 @@ 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(1.4826*scipy.stats.median_abs_deviation(deltaMag[goodDelta]),
+                          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 - zp) <= self.config.magnitudeOutlierRejectionNSigma*zpSigma)]
+
+        nOutlier = nMatch - goodStars.size
+        self.log.info("Removing %d of %d magnitude outliers.", nOutlier, nMatch)
 
         matchesOut = []
-        for matchInd in good:
+        for matchInd in goodStars:
             matchesOut.append(matchesIn[matchInd])
 
         return matchesOut

tests/test_astrometryTask.py

@@ -170,7 +170,7 @@ def doTest(self, pixelsToTanPixels, order=3):
             sourceCat=sourceCat,
             exposure=self.exposure,
         )
-        self.assertEqual(len(resultsMagOutlierRejection.matches), len(resultsRefSelect.matches))
+        self.assertLess(len(resultsMagOutlierRejection.matches), len(resultsRefSelect.matches))
         config.doMagnitudeOutlierRejection = False
 
         # try again, but without fitting the WCS, no reference selector
@@ -222,6 +222,18 @@ def makeSourceCat(self, distortedWcs):
             src.set(sourceCentroidKey, refObj.get(refCentroidKey))
             src.set(sourceInstFluxKey, refObj.get(refFluxRKey))
             src.set(sourceInstFluxErrKey, refObj.get(refFluxRKey)/100)
+
+        # Make a few outlier fluxes for the outlier rejection tests.
+        # Otherwise, these fluxes are not used.
+        unResolved, = np.where(~refCat['resolved'])
+
+        sourceCat[sourceInstFluxKey][unResolved[0]] = 0.0
+        sourceCat[sourceInstFluxKey][unResolved[1]] = -1.0
+        sourceCat[sourceInstFluxKey][unResolved[2]] *= 1000.0
+        refCat[refFluxRKey][unResolved[3]] = 0.0
+        refCat[refFluxRKey][unResolved[4]] = -1.0
+        refCat[refFluxRKey][unResolved[5]] *= 1000.0
+
         return sourceCat