Improved AstrometryTask iteration and FitTanWcsTask fitting

AstrometryTask now iterates the match and fit WCS cycle in a smarter way:
- It checks if each iteration has improved the fit and stops if not
  or if the fit is so good there is no point continuing.
- If a WCS fit fails it either gives up (if it's the first fit)
  or uses the previous iteration
- It reduces the maximum matching distance for each iteration
  based on the measured quality of the fit.
Improved FitTanSipWcsTask in two ways:
- It iterates the fit (since the existing C++ code fits X first, then Y,
  which is something we plan to fix at some point).
- If the fit is terrible it raises an exception.

python/lsst/meas/astrom/astrometry.py

@@ -1,10 +1,9 @@
 from __future__ import absolute_import, division, print_function
 
 from lsst.daf.base import PropertyList
-from lsst.afw.image import ExposureF
 from lsst.afw.image.utils import getDistortedWcs
-from lsst.afw.table import Point2DKey
-from lsst.afw.geom import Box2D
+from lsst.afw.geom import Box2D, radToArcsec
+import lsst.afw.math as afwMath
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
 from .loadAstrometryNetObjects import LoadAstrometryNetObjectsTask
@@ -35,7 +34,7 @@ class AstrometryConfig(pexConfig.Config):
         doc = "maximum number of iterations of match sources and fit WCS; " +
             "ignored if forceKnownWcs True",
         dtype = int,
-        default = 3,
+        default = 6,
         min = 1,
     )
 
@@ -196,40 +195,53 @@ def solve(self, sourceCat, bbox, initWcs, filterName=None, calib=None, exposure=
 
         res = None
         wcs = initWcs
+        maxMatchDistArcSec = None
         for i in range(self.config.maxIter):
-            tryRes = self._matchAndFitWcs( # refCat, sourceCat, refFluxField, bbox, wcs, exposure=None
-                refCat = loadRes.refCat,
-                sourceCat = sourceCat,
-                refFluxField = loadRes.fluxField,
-                bbox = bbox,
-                wcs = wcs,
-                exposure = exposure,
-            )
+            try:
+                tryRes = self._matchAndFitWcs( # refCat, sourceCat, refFluxField, bbox, wcs, exposure=None
+                    refCat = loadRes.refCat,
+                    sourceCat = sourceCat,
+                    refFluxField = loadRes.fluxField,
+                    bbox = bbox,
+                    wcs = wcs,
+                    exposure = exposure,
+                    maxMatchDistArcSec = maxMatchDistArcSec,
+                )
+            except Exception as e:
+                # if we have had a succeessful iteration then use that; otherwise fail
+                if i > 0:
+                    self.log.info("Fit WCS iter %d failed; using previous iteration: %s" % (i, e))
+                    break
+                else:
+                    raise
 
             if self.config.forceKnownWcs:
                 # run just once; note that the number of matches has already logged
                 res = tryRes
                 break
 
-            self.log.logdebug(
-                "Fit WCS iter %s: %s matches; median scatter = %g arcsec" % \
-                    (i, len(tryRes.matches), tryRes.scatterOnSky.asArcseconds()),
-            )
+            distRadList = [match.distance for match in tryRes.matches]
+            distRadStats = afwMath.makeStatistics(distRadList, afwMath.MEANCLIP | afwMath.STDEVCLIP)
+            distArcsecMean = radToArcsec(distRadStats.getValue(afwMath.MEANCLIP))
+            distArcsecStdDev = radToArcsec(distRadStats.getValue(afwMath.STDEVCLIP))
 
-            if res is not None and not self.config.forceKnownWcs:
-                if len(tryRes.matches) < len(res.matches):
-                    self.log.info(
-                        "Fit WCS: use iter %s because it had more matches than the next iter: %s vs. %s" % \
-                        (i-1, len(res.matches), len(tryRes.matches)))
-                    break
-                if len(tryRes.matches) == len(res.matches) and tryRes.scatterOnSky >= res.scatterOnSky:
-                    self.log.info(
-            "Fit WCS: use iter %s because it had less scatter than the next iter: %g vs. %g arcsec" % \
-                        (i-1, res.scatterOnSky.asArcseconds(), tryRes.scatterOnSky.asArcseconds()))
+            self.log.info("Fit WCS iter %d: %d matches; scatter = %0.3f +- %0.3f arcsec" %
+                    (i, len(tryRes.matches), distArcsecMean, distArcsecStdDev))
+
+            newMaxMatchDistArcSec = distArcsecMean + 2.0*distArcsecStdDev
+            if maxMatchDistArcSec is not None:
+                if newMaxMatchDistArcSec >= maxMatchDistArcSec:
+                    self.log.warn(
+                        "Iteration %d had no better mean + 2sigma scatter; using previous iteration" % (i,))
                     break
 
+            maxMatchDistArcSec = newMaxMatchDistArcSec
             res = tryRes
             wcs = res.wcs
+            if newMaxMatchDistArcSec < 0.001:
+                self.log.info(
+                    "Iteration %d had mean + 2sigma scatter < 0.001 arcsec; that's good enough" % (i,))
+                break
 
         return pipeBase.Struct(
             refCat = loadRes.refCat,
@@ -241,13 +253,16 @@ def solve(self, sourceCat, bbox, initWcs, filterName=None, calib=None, exposure=
         )
 
     @pipeBase.timeMethod
-    def _matchAndFitWcs(self, refCat, sourceCat, refFluxField, bbox, wcs, exposure=None):
+    def _matchAndFitWcs(self, refCat, sourceCat, refFluxField, bbox, wcs, maxMatchDistArcSec=None,
+        exposure=None):
         """!Match sources to reference objects and fit a WCS
 
         @param[in] refCat  catalog of reference objects
         @param[in] sourceCat  catalog of sourceCat detected on the exposure (an lsst.afw.table.SourceCatalog)
         @param[in] bbox  bounding box of exposure (an lsst.afw.geom.Box2I)
         @param[in] wcs  initial guess for WCS of exposure (an lsst.afw.image.Wcs)
+        @param[in] maxMatchDistArcSec  maximum distance between reference objects and sources (arcsec);
+            if None then use the matcher's default
         @param[in] exposure  exposure whose WCS is to be fit, or None; used only for the debug display
 
         @return an lsst.pipe.base.Struct with these fields:
@@ -263,7 +278,9 @@ def _matchAndFitWcs(self, refCat, sourceCat, refFluxField, bbox, wcs, exposure=N
             sourceCat = sourceCat,
             wcs = wcs,
             refFluxField = refFluxField,
+            maxMatchDistArcSec = maxMatchDistArcSec,
         )
+        self.log.info("Found %s matches" % (len(matchRes.matches),))
         if debug.display:
             frame = int(debug.frame)
             displayAstrometry(
@@ -288,7 +305,7 @@ def _matchAndFitWcs(self, refCat, sourceCat, refFluxField, bbox, wcs, exposure=N
             fitWcs = fitRes.wcs
             scatterOnSky = fitRes.scatterOnSky
         else:
-            self.log.info("Not fitting WCS (forceKnownWcs true); %s matches" % (len(matchRes.matches),))
+            self.log.info("Not fitting WCS (forceKnownWcs true)")
             fitWcs = wcs
             scatterOnSky = None
         if debug.display:

python/lsst/meas/astrom/fitTanSipWcs.py

@@ -11,9 +11,23 @@
 
 class FitTanSipWcsConfig(pexConfig.Config):
     order = pexConfig.RangeField(
-        doc = "order of SIP polynomial (0 for pure TAN WCS)",
+        doc = "order of SIP polynomial",
+        dtype = int,
+        default = 4,
+        min = 0,
+    )
+    numIter = pexConfig.RangeField(
+        doc = "number of iterations of fitter (which fits X and Y separately, and so benefits from " + \
+            "a few iterations",
         dtype = int,
         default = 3,
+        min = 1,
+    )
+    maxScatterArcsec = pexConfig.RangeField(
+        doc = "maximum median scatter of a WCS fit beyond which the fit fails (arcsec); " +
+            "be generous, as this is only intended to catch catastrophic failures",
+        dtype = float,
+        default = 10,
         min = 0,
     )
 
@@ -101,8 +115,10 @@ def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None):
         """
         if bbox is None:
             bbox = afwGeom.Box2I()
-        sipObject = makeCreateWcsWithSip(matches, initWcs, self.config.order, bbox)
-        wcs = sipObject.getNewWcs()
+        wcs = initWcs
+        for i in range(self.config.numIter):
+            sipObject = makeCreateWcsWithSip(matches, wcs, self.config.order, bbox)
+            wcs = sipObject.getNewWcs()
 
         if refCat is not None:
             self.log.info("Updating centroids in refCat")
@@ -121,9 +137,16 @@ def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None):
         self.log.info("Updating distance in match list")
         setMatchDistance(matches)
 
+        scatterOnSky = sipObject.getScatterOnSky()
+
+        if scatterOnSky.asArcseconds() > self.config.maxScatterArcsec:
+            raise pipeBase.TaskError(
+                "Fit failed: median scatter on sky = %0.3f arcsec > %0.3f config.maxScatterArcsec" %
+                (scatterOnSky.asArcseconds(), self.config.maxScatterArcsec))
+
         return pipeBase.Struct(
             wcs = wcs,
-            scatterOnSky = sipObject.getScatterOnSky(),
+            scatterOnSky = scatterOnSky,
         )
 
     @staticmethod

tests/testAstrometryTask.py

@@ -26,7 +26,7 @@
 import os
 import unittest
 
-import numpy
+import numpy as np
 
 import eups
 import lsst.utils.tests as utilsTests
@@ -116,7 +116,7 @@ def doTest(self, pixelsToTanPixels, order=3):
 
             angSep = refCoord.angularSeparation(srcCoord)
             maxAngSep = max(maxAngSep, angSep)
-            self.assertLess(refCoord.angularSeparation(srcCoord), 0.0025 * afwGeom.arcseconds)
+            self.assertLess(refCoord.angularSeparation(srcCoord).asArcseconds(), 0.0025)
 
             pixSep = math.hypot(*(srcPixPos-refPixPos))
             maxPixSep = max(maxPixSep, pixSep)
@@ -136,8 +136,11 @@ def doTest(self, pixelsToTanPixels, order=3):
         self.assertTrue(resultsNoFit.wcs is distortedWcs)
         self.assertTrue(resultsNoFit.initWcs is distortedWcs)
         self.assertTrue(resultsNoFit.scatterOnSky is None)
-        # fitting may find a few more matches, since it matches again after fitting the WCS
-        self.assertTrue(0 <= len(results.matches) - len(resultsNoFit.matches) < len(results.matches) * 0.1)
+
+        # fitting should improve the quality of the matches
+        meanFitDist = np.mean([match.distance for match in results.matches])
+        meanNoFitDist = np.mean([match.distance for match in results.matches])
+        self.assertLessEqual(meanFitDist, meanNoFitDist)
 
     def makeSourceCat(self, distortedWcs):
         """Make a source catalog by reading the position reference stars and distorting the positions
@@ -147,17 +150,18 @@ def makeSourceCat(self, distortedWcs):
         loadRes = loader.loadPixelBox(bbox=self.bbox, wcs=distortedWcs, filterName="r")
         refCat = loadRes.refCat
         refCentroidKey = afwTable.Point2DKey(refCat.schema["centroid"])
+        refFluxRKey = refCat.schema["r_flux"].asKey()
 
         sourceSchema = afwTable.SourceTable.makeMinimalSchema()
         measBase.SingleFrameMeasurementTask(schema=sourceSchema) # expand the schema
         sourceCat = afwTable.SourceCatalog(sourceSchema)
         sourceCentroidKey = afwTable.Point2DKey(sourceSchema["slot_Centroid"])
+        sourceFluxKey = sourceSchema["slot_ApFlux_flux"].asKey()
 
         for refObj in refCat:
             src = sourceCat.addNew()
-            refPos = refObj.get(refCentroidKey)
-            srcPos = refPos
-            src.set(sourceCentroidKey, srcPos)
+            src.set(sourceCentroidKey, refObj.get(refCentroidKey))
+            src.set(sourceFluxKey, refObj.get(refFluxRKey))
         return sourceCat
 
     def assertWcssAlmostEqual(self, wcs0, wcs1, bbox, maxSkyErr=0.01 * afwGeom.arcseconds, maxPixErr = 0.02,
@@ -177,8 +181,8 @@ def assertWcssAlmostEqual(self, wcs0, wcs1, bbox, maxSkyErr=0.01 * afwGeom.arcse
         @throw AssertionError if the two WCSs do not match sufficiently closely
         """
         bboxd = afwGeom.Box2D(bbox)
-        xList = numpy.linspace(bboxd.getMinX(), bboxd.getMaxX(), nx)
-        yList = numpy.linspace(bboxd.getMinY(), bboxd.getMaxY(), ny)
+        xList = np.linspace(bboxd.getMinX(), bboxd.getMaxX(), nx)
+        yList = np.linspace(bboxd.getMinY(), bboxd.getMaxY(), ny)
         maxSkyErrPixPos = [afwGeom.Angle(0), None]
         maxPixErrSkyPos = [0, None]
         for x in xList:

tests/testFitTanSipWcsTask.py

@@ -55,7 +55,7 @@ class BaseTestCase(unittest.TestCase):
 
     def setUp(self):
         crval = afwCoord.IcrsCoord(afwGeom.PointD(44., 45.))
-        crpix = afwGeom.PointD(0, 0)
+        crpix = afwGeom.Point2D(15000, 4000)
 
         arcsecPerPixel = 1/3600.0
         CD11 = arcsecPerPixel
@@ -66,7 +66,7 @@ def setUp(self):
         self.tanWcs = afwImage.makeWcs(crval, crpix, CD11, CD12, CD21, CD22)
         self.loadData()
 
-    def loadData(self, rangePix=3000, numPoints=5):
+    def loadData(self, rangePix=3000, numPoints=25):
         """Load catalogs and make the match list
 
         This is a separate function so data can be reloaded if fitting more than once
@@ -117,35 +117,41 @@ def tearDown(self):
 
     def testTrivial(self):
         """Add no distortion"""
-        self.doTest("testTrivial", lambda x, y: (x, y))
+        for order in (2, 4, 6):
+            self.doTest("testTrivial", lambda x, y: (x, y), order=order)
 
     def testOffset(self):
         """Add an offset"""
-        self.doTest("testOffset", lambda x, y: (x + 5, y + 7))
+        for order in (2, 4, 6):
+            self.doTest("testOffset", lambda x, y: (x + 5, y + 7), order=order)
 
     def testLinearX(self):
         """Scale x, offset y"""
-        self.doTest("testLinearX", lambda x, y: (2*x, y + 7))
+        for order in (2, 6):
+            self.doTest("testLinearX", lambda x, y: (2*x, y + 7), order=order)
 
     def testLinearXY(self):
         """Scale x and y"""
         self.doTest("testLinearXY", lambda x, y: (2*x, 3*y))
 
     def testLinearYX(self):
         """Add an offset to each point; scale in y and x"""
-        self.doTest("testLinearYX", lambda x, y: (x + 0.2*y, y + 0.3*x))
+        for order in (2, 6):
+            self.doTest("testLinearYX", lambda x, y: (x + 0.2*y, y + 0.3*x), order=order)
 
     def testQuadraticX(self):
         """Add quadratic distortion in x"""
-        self.doTest("testQuadraticX", lambda x, y: (x + 1e-5*x**2, y), order=4, specifyBBox=True)
+        for order in (4, 5):
+            self.doTest("testQuadraticX", lambda x, y: (x + 1e-5*x**2, y), order=order)
 
     def testRadial(self):
         """Add radial distortion"""
         radialTransform = afwGeom.RadialXYTransform([0, 1.01, 1e-8])
         def radialDistortion(x, y):
             x, y = radialTransform.forwardTransform(afwGeom.Point2D(x, y))
             return (x, y)
-        self.doTest("testRadial", radialDistortion)
+        for order in (4, 5, 6):
+            self.doTest("testRadial", radialDistortion, order=order)
 
     def checkResults(self, fitRes, catsUpdated):
         """Check results
@@ -179,17 +185,22 @@ def checkResults(self, fitRes, catsUpdated):
             distErr = abs(dist - angSep)
             maxDistErr = max(maxDistErr, distErr)
             maxAngSep = max(maxAngSep, angSep)
-            self.assertLess(angSep, 0.001 * afwGeom.arcseconds)
+            self.assertLess(angSep.asArcseconds(), 0.001)
 
             pixSep = math.hypot(*(srcPixPos - refPixPos))
             maxPixSep = max(maxPixSep, pixSep)
             self.assertLess(pixSep, 0.001)
 
         print("max angular separation = %0.4f arcsec" % (maxAngSep.asArcseconds(),))
         print("max pixel separation = %0.3f" % (maxPixSep,))
-        self.assertLess(maxDistErr.asArcseconds(), 1e-7)
+        if catsUpdated:
+            allowedDistErr = 1e-7
+        else:
+            allowedDistErr = 0.001
+        self.assertLess(maxDistErr.asArcseconds(), allowedDistErr,
+            "Computed distance in match list is off by %s arcsec" % (maxDistErr.asArcseconds(),))
 
-    def doTest(self, name, func, order=3, specifyBBox=False, doPlot=False):
+    def doTest(self, name, func, order=3, numIter=4, specifyBBox=False, doPlot=False):
         """Apply func(x, y) to each source in self.sourceCat, then fit and check the resulting WCS
         """
         bbox = afwGeom.Box2I()
@@ -200,11 +211,13 @@ def doTest(self, name, func, order=3, specifyBBox=False, doPlot=False):
             src.set(self.srcCentroidKey, distortedPos)
             bbox.include(afwGeom.Point2I(afwGeom.Point2I(distortedPos)))
 
-        if specifyBBox:
-            sipObject = makeCreateWcsWithSip(self.matches, self.tanWcs, order, bbox)
-        else:
-            sipObject = makeCreateWcsWithSip(self.matches, self.tanWcs, order)
-        tanSipWcs = sipObject.getNewWcs()
+        tanSipWcs = self.tanWcs
+        for i in range(numIter):
+            if specifyBBox:
+                sipObject = makeCreateWcsWithSip(self.matches, tanSipWcs, order, bbox)
+            else:
+                sipObject = makeCreateWcsWithSip(self.matches, tanSipWcs, order)
+            tanSipWcs = sipObject.getNewWcs()
         setMatchDistance(self.matches)
         fitRes = lsst.pipe.base.Struct(
             wcs = tanSipWcs,
@@ -223,6 +236,7 @@ def doTest(self, name, func, order=3, specifyBBox=False, doPlot=False):
 
             fitterConfig = FitTanSipWcsTask.ConfigClass()
             fitterConfig.order = order
+            fitterConfig.numIter = numIter
             fitter = FitTanSipWcsTask(config=fitterConfig)
             self.loadData()
             if specifyBBox: