Set distance field in matches in matcher and WCS fitter

python/lsst/meas/astrom/__init__.py

@@ -28,6 +28,7 @@
 from .anetAstrometry import *
 from .loadAstrometryNetObjects import *
 from .matchOptimisticB import *
+from .setMatchDistance import *
 
 from . import catalogStarSelector
 

python/lsst/meas/astrom/fitTanSipWcs.py

@@ -4,6 +4,7 @@
 import lsst.afw.table as afwTable
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
+from .setMatchDistance import setMatchDistance
 from .sip import makeCreateWcsWithSip
 
 __all__ = ["FitTanSipWcsTask", "FitTanSipWcsConfig"]
@@ -74,19 +75,23 @@ def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None):
 
         @param[in,out] matches  a list of reference object/source matches
             (an lsst::afw::table::ReferenceMatchVector)
-            The centroids of the reference objects in matches may be updated for the new WCS,
-            but it is strongly recommended that you provide the refCat and sourceCat arguments
-            so that all matches are fully updated, as well as any sources or reference objects
-            not in matches.
+            The following fields are read:
+            - match.first (reference object) coord
+            - match.second (source) centroid
+            The following fields are written:
+            - match.first (reference object) centroid,
+            - match.second (source) centroid
+            - match.distance (on sky separation, in radians)
         @param[in] initWcs  initial WCS
         @param[in] bbox  the region over which the WCS will be valid (an lsst:afw::geom::Box2I);
             if None or an empty box then computed from matches
         @param[in,out] refCat  reference object catalog, or None.
             If provided then all centroids are updated with the new WCS,
-            otherwise the centroids in matches might be updated but the others will not be touched.
+            otherwise only the centroids for ref objects in matches are updated.
             Required fields are "centroid_x", "centroid_y" and "coord".
         @param[in,out] sourceCat  source catalog, or None.
-            If provided then coords are updated with the new WCS.
+            If provided then coords are updated with the new WCS;
+            otherwise only the coords for sources in matches are updated.
             Required fields are "slot_Centroid_x", "slot_Centroid_y" and "coord".
 
         @return an lsst.pipe.base.Struct with the following fields:
@@ -100,30 +105,47 @@ def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None):
         wcs = sipObject.getNewWcs()
 
         if refCat is not None:
-            self.updateRefCat(wcs, refCat)
+            self.log.info("Updating centroids in refCat")
+            self.updateRefCentroids(wcs, refList=refCat)
+        else:
+            self.log.warning("Updating reference object centroids in match list; refCat is None")
+            self.updateRefCentroids(wcs, refList=[match.first for match in matches])
 
         if sourceCat is not None:
-            self.updateSrcCat(wcs, sourceCat)
+            self.log.info("Updating coords in sourceCat")
+            self.updateSourceCoords(wcs, sourceList=sourceCat)
+        else:
+            self.log.warning("Updating source coords in match list; sourceCat is None")
+            self.updateSourceCoords(wcs, sourceList=[match.second for match in matches])
+
+        self.log.info("Updating distance in match list")
+        setMatchDistance(matches)
 
         return pipeBase.Struct(
             wcs = wcs,
             scatterOnSky = sipObject.getScatterOnSky(),
         )
 
     @staticmethod
-    def updateRefCat(wcs, refCat):
-        """Update centroids in a reference catalog, given a WCS
+    def updateRefCentroids(wcs, refList):
+        """Update centroids in a collection of reference objects, given a WCS
         """
-        coordKey = refCat.schema["coord"].asKey()
-        centroidKey = afwTable.Point2DKey(refCat.schema["centroid"])
-        for refObj in refCat:
+        if len(refList) < 1:
+            return
+        schema = refList[0].schema
+        coordKey = schema["coord"].asKey()
+        centroidKey = afwTable.Point2DKey(schema["centroid"])
+        for refObj in refList:
             refObj.set(centroidKey, wcs.skyToPixel(refObj.get(coordKey)))
 
     @staticmethod
-    def updateSrcCat(wcs, sourceCat):
-        """Update coords in a source catalog, given a WCS
+    def updateSourceCoords(wcs, sourceList):
+        """Update coords in a collection of sources, given a WCS
         """
-        srcCoordKey = sourceCat.schema["coord"].asKey()
-        for src in sourceCat:
+        if len(sourceList) < 1:
+            return
+        schema = sourceList[1].schema
+        srcCoordKey = schema["coord"].asKey()
+        for src in sourceList:
             src.set(srcCoordKey, wcs.pixelToSky(src.getCentroid()))
 

python/lsst/meas/astrom/matchOptimisticB.py

@@ -5,6 +5,7 @@
 import lsst.afw.table as afwTable
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
+from .setMatchDistance import setMatchDistance
 from .astromLib import matchOptimisticB, MatchOptimisticBControl
 
 __all__ = ["MatchOptimisticBTask", "MatchOptimisticBConfig"]
@@ -323,6 +324,8 @@ def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numCleanSources, minMat
         @param[in] minMatchedPairs  minimum number of matches
         @param[in] sourceInfo  SourceInfo for the sourceCat
         @param[in] verbose  true to print diagnostic information to std::cout
+
+        @return a list of matches, an instance of lsst.afw.table.ReferenceMatch
         """
         numSources = len(sourceCat)
         posRefBegInd = numCleanSources - numSources
@@ -352,4 +355,6 @@ def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numCleanSources, minMat
                         verbose,
                     )
                     if matches is not None and len(matches) != 0:
+                        setMatchDistance(matches)
                         return matches
+

python/lsst/meas/astrom/setMatchDistance.py

@@ -0,0 +1,42 @@
+from __future__ import absolute_import, division, print_function
+# 
+# LSST Data Management System
+# Copyright 2008, 2009, 2010 LSST Corporation.
+# 
+# This product includes software developed by the
+# LSST Project (http://www.lsst.org/).
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+# 
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the LSST License Statement and 
+# the GNU General Public License along with this program.  If not, 
+# see <http://www.lsstcorp.org/LegalNotices/>.
+#
+__all__ = ["setMatchDistance"]
+
+def setMatchDistance(matches):
+    """Set the distance field of the matches in a match list to the distance in radians on the sky
+
+    @warning the coord field of the source in each match must be correct
+
+    @param[in,out] matches  a list of matches, an instance of lsst.afw.table.ReferenceMatch
+        reads the coord field of the source and reference object of each match
+        writes the distance field of each match
+    """
+    if len(matches) < 1:
+        return
+
+    sourceCoordKey = matches[0].first.schema["coord"].asKey()
+    refObjCoordKey = matches[0].second.schema["coord"].asKey()
+    for match in matches:
+        sourceCoord = match.first.get(sourceCoordKey)
+        refObjCoord = match.second.get(refObjCoordKey)
+        match.distance = refObjCoord.angularSeparation(sourceCoord).asRadians()

tests/testFitTanSipWcsTask.py

@@ -39,7 +39,7 @@
 import lsst.afw.table as afwTable
 from lsst.meas.algorithms import LoadReferenceObjectsTask
 from lsst.meas.base import SingleFrameMeasurementTask
-from lsst.meas.astrom import FitTanSipWcsTask
+from lsst.meas.astrom import FitTanSipWcsTask, setMatchDistance
 from lsst.meas.astrom.sip import makeCreateWcsWithSip
 
 class BaseTestCase(unittest.TestCase):
@@ -64,10 +64,14 @@ def setUp(self):
         CD22 = arcsecPerPixel
 
         self.tanWcs = afwImage.makeWcs(crval, crpix, CD11, CD12, CD21, CD22)
+        self.loadData()
 
-        S = 3000
-        N = 5
+    def loadData(self, rangePix=3000, numPoints=5):
+        """Load catalogs and make the match list
 
+        This is a separate function so data can be reloaded if fitting more than once
+        (each time a WCS is fit it may update the source catalog, reference catalog and match list)
+        """
         if self.MatchClass == afwTable.ReferenceMatch:
             refSchema = LoadReferenceObjectsTask.makeMinimalSchema(
                 filterNameList = ["r"], addFluxSigma=True, addIsPhotometric=True)
@@ -84,22 +88,17 @@ def setUp(self):
         self.srcCentroidKey_xSigma = srcSchema["slot_Centroid_xSigma"].asKey()
         self.srcCentroidKey_ySigma = srcSchema["slot_Centroid_ySigma"].asKey()
         self.sourceCat = afwTable.SourceCatalog(srcSchema)
-        self.origSourceCat = afwTable.SourceCatalog(srcSchema) # undistorted copy
-        self.matches = []
 
+        self.matches = []
 
-        for i in numpy.linspace(0., S, N):
-            for j in numpy.linspace(0., S, N):
+        for i in numpy.linspace(0., rangePix, numPoints):
+            for j in numpy.linspace(0., rangePix, numPoints):
                 src = self.sourceCat.addNew()
-                origSrc = self.origSourceCat.addNew()
                 refObj = self.refCat.addNew()
 
                 src.set(self.srcCentroidKey, afwGeom.Point2D(i, j))
                 src.set(self.srcCentroidKey_xSigma, 0.1)
                 src.set(self.srcCentroidKey_ySigma, 0.1)
-                origSrc.set(self.srcCentroidKey, afwGeom.Point2D(i, j))
-                origSrc.set(self.srcCentroidKey_xSigma, 0.1)
-                origSrc.set(self.srcCentroidKey_ySigma, 0.1)
 
                 c = self.tanWcs.pixelToSky(afwGeom.Point2D(i, j))
                 refObj.setCoord(c)
@@ -112,7 +111,6 @@ def setUp(self):
 
     def tearDown(self):
         del self.refCat
-        del self.origSourceCat
         del self.sourceCat
         del self.matches
         del self.tanWcs
@@ -165,7 +163,8 @@ def checkResults(self, fitRes, catsUpdated):
         refCoordKey = self.refCat.schema["coord"].asKey()
         if catsUpdated:
             refCentroidKey = afwTable.Point2DKey(self.refCat.schema["centroid"])
-        for refObj, src, d in self.matches:
+        maxDistErr = afwGeom.Angle(0)
+        for refObj, src, distRad in self.matches:
             srcPixPos = src.get(self.srcCentroidKey)
             refCoord = refObj.get(refCoordKey)
             if catsUpdated:
@@ -176,6 +175,9 @@ def checkResults(self, fitRes, catsUpdated):
                 srcCoord = tanSipWcs.pixelToSky(srcPixPos)
 
             angSep = refCoord.angularSeparation(srcCoord)
+            dist = distRad*afwGeom.radians
+            distErr = abs(dist - angSep)
+            maxDistErr = max(maxDistErr, distErr)
             maxAngSep = max(maxAngSep, angSep)
             self.assertLess(angSep, 0.001 * afwGeom.arcseconds)
 
@@ -185,6 +187,7 @@ def checkResults(self, fitRes, catsUpdated):
 
         print("max angular separation = %0.4f arcsec" % (maxAngSep.asArcseconds(),))
         print("max pixel separation = %0.3f" % (maxPixSep,))
+        self.assertLess(maxDistErr.asArcseconds(), 1e-7)
 
     def doTest(self, name, func, order=3, specifyBBox=False, doPlot=False):
         """Apply func(x, y) to each source in self.sourceCat, then fit and check the resulting WCS
@@ -202,6 +205,7 @@ def doTest(self, name, func, order=3, specifyBBox=False, doPlot=False):
         else:
             sipObject = makeCreateWcsWithSip(self.matches, self.tanWcs, order)
         tanSipWcs = sipObject.getNewWcs()
+        setMatchDistance(self.matches)
         fitRes = lsst.pipe.base.Struct(
             wcs = tanSipWcs,
             scatterOnSky = sipObject.getScatterOnSky(),
@@ -213,9 +217,14 @@ def doTest(self, name, func, order=3, specifyBBox=False, doPlot=False):
         self.checkResults(fitRes, catsUpdated=False)
 
         if self.MatchClass == afwTable.ReferenceMatch:
+            # reset source coord and reference centroid based on initial WCS
+            FitTanSipWcsTask.updateRefCentroids(wcs=self.tanWcs, refList = self.refCat)
+            FitTanSipWcsTask.updateSourceCoords(wcs=self.tanWcs, sourceList = self.sourceCat)
+
             fitterConfig = FitTanSipWcsTask.ConfigClass()
             fitterConfig.order = order
             fitter = FitTanSipWcsTask(config=fitterConfig)
+            self.loadData()
             if specifyBBox:
                 fitRes = fitter.fitWcs(
                     matches = self.matches,

tests/testMatchOptimisticB.py

@@ -102,11 +102,19 @@ def singleTestInstance(self, filename, distortFunc, doPlot=False):
         if doPlot:
             self.plotMatches(matches)
         self.assertEqual(len(matches), 183)
+
+        refCoordKey = refCat.schema["coord"].asKey()
+        srcCoordKey = sourceCat.schema["coord"].asKey()
         numErrors = 0
-        for match in matches:
-            if match.first.getId() != match.second.getId():
-                refObj = match.first
-                source = match.second
+        maxDistErr = afwGeom.Angle(0)
+        for refObj, source, distRad in matches:
+            sourceCoord = source.get(srcCoordKey)
+            refCoord = refObj.get(refCoordKey)
+            predDist = sourceCoord.angularSeparation(refCoord)
+            distErr = abs(predDist - distRad*afwGeom.radians)
+            maxDistErr = max(distErr, maxDistErr)
+
+            if refObj.getId() != source.getId():
                 refCentroid = refObj.get("centroid")
                 sourceCentroid = source.getCentroid()
                 radius = math.hypot(*(refCentroid - sourceCentroid))
@@ -115,6 +123,7 @@ def singleTestInstance(self, filename, distortFunc, doPlot=False):
                     source.getId(), sourceCentroid, radius))
         print("num match errors=", numErrors)
         self.assertLess(numErrors, 3)
+        self.assertLess(maxDistErr.asArcseconds(), 1e-7)
 
     def computePosRefCatalog(self, sourceCat):
         """Generate a position reference catalog from a source catalog