Merge pull request #11 from lsst/tickets/DM-15008

DM-15008: anetAstrometry.py uses self.distortionContext, which does not exist

.gitignore

@@ -9,6 +9,7 @@ _build.*
 .cache
 .pytest_cache
 pytest_session.txt
+.coverage
 config.log
 doc/html
 doc/doxygen.conf

python/lsst/meas/extensions/astrometryNet/anetAstrometry.py

@@ -33,7 +33,7 @@
 import lsst.pex.exceptions
 import lsst.afw.geom as afwGeom
 from lsst.afw.cameraGeom import PIXELS, TAN_PIXELS
-from lsst.afw.table import Point2DKey, CovarianceMatrix2fKey
+from lsst.afw.table import Point2DKey, CovarianceMatrix2fKey, updateSourceCoords
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
 from lsst.meas.astrom import displayAstrometry
@@ -140,6 +140,7 @@ def DebugInfo(name):
     into your debug.py file and run photoCalTask.py with the \c --debug flag.
     """
     ConfigClass = ANetAstrometryConfig
+    _DefaultName = "astrometricSolver"
 
     def __init__(self, schema, refObjLoader=None, **kwds):
         """!Create the astrometric calibration task.  Most arguments are simply passed onto pipe.base.Task.
@@ -224,8 +225,7 @@ def solve(self, exposure, sourceCat):
 
         \note ignores config.forceKnownWcs
         """
-        with self.distortionContext(sourceCat=sourceCat, exposure=exposure) as bbox:
-            results = self._astrometry(sourceCat=sourceCat, exposure=exposure, bbox=bbox)
+        results = self._astrometry(sourceCat=sourceCat, exposure=exposure)
 
         if results.matches:
             self.refitWcs(sourceCat=sourceCat, exposure=exposure, matches=results.matches)
@@ -320,7 +320,7 @@ def loadAndMatch(self, exposure, sourceCat, bbox=None):
         matchMeta = astrom.getMatchMetadata()
         if matches is None or len(matches) == 0:
             raise RuntimeError("No astrometric matches")
-        self.log.info("%d astrometric matches" % (len(matches)))
+        self.log.info("%d astrometric matches", len(matches))
 
         if self._display:
             frame = lsstDebug.Info(__name__).frame
@@ -362,7 +362,7 @@ def _astrometry(self, sourceCat, exposure, bbox=None):
         matchMeta = astrom.getMatchMetadata()
         if matches is None or len(matches) == 0:
             raise RuntimeError("No astrometric matches")
-        self.log.info("%d astrometric matches" % (len(matches)))
+        self.log.info("%d astrometric matches", len(matches))
 
         # Note that this is the Wcs for the provided positions, which may be distorted
         exposure.setWcs(astrom.getWcs())
@@ -434,17 +434,14 @@ def fitWcs(initialWcs, title=None):
                 wcs, scatter = fitWcs(wcs, title="Final astrometry")
 
             except lsst.pex.exceptions.LengthError as e:
-                self.log.warn("Unable to fit SIP: %s" % e)
+                self.log.warn("Unable to fit SIP: %s", e)
 
-            self.log.info("Astrometric scatter: %f arcsec (%s non-linear terms, %d matches, %d rejected)" %
-                          (scatter.asArcseconds(), "with" if wcs.hasDistortion() else "without",
-                           len(matches), numRejected))
+            self.log.info("Astrometric scatter: %f arcsec (%d matches, %d rejected)",
+                          scatter.asArcseconds(), len(matches), numRejected)
             exposure.setWcs(wcs)
 
             # Apply WCS to sources
-            for index, source in enumerate(sourceCat):
-                sky = wcs.pixelToSky(source.getX(), source.getY())
-                source.setCoord(sky)
+            updateSourceCoords(wcs, sourceCat)
         else:
             self.log.warn("Not calculating a SIP solution; matches may be suspect")
 

python/lsst/meas/extensions/astrometryNet/anetBasicAstrometry.py

@@ -46,6 +46,7 @@
 from .loadAstrometryNetObjects import LoadAstrometryNetObjectsTask, LoadMultiIndexes
 from lsst.meas.astrom import displayAstrometry, makeMatchStatisticsInRadians
 import lsst.meas.astrom.sip as astromSip
+from . import cleanBadPoints
 
 
 class InitialAstrometry(object):
@@ -170,7 +171,7 @@ class ANetBasicAstrometryConfig(LoadAstrometryNetObjectsTask.ConfigClass):
         min=0.0,
     )
     cleaningParameter = RangeField(
-        doc="Sigma-clipping parameter in sip/cleanBadPoints.py",
+        doc="Sigma-clipping parameter in cleanBadPoints.py",
         dtype=float,
         default=3.0,
         min=0.0,
@@ -844,7 +845,7 @@ def _getMatchList(self, sourceCat, refCat, wcs):
                 self.loginfo('Reference range: RA [%.3f, %.3f], Dec [%.3f, %.3f]' %
                              (min(R2), max(R2), min(D2), max(D2)))
             raise RuntimeError('No matches found between image and catalogue')
-        matches = astromSip.cleanBadPoints.clean(matches, wcs, nsigma=clean)
+        matches = cleanBadPoints.clean(matches, wcs, nsigma=clean)
         return matches
 
     def getColumnName(self, filterName, columnMap, default=None):

python/lsst/meas/extensions/astrometryNet/cleanBadPoints.py

@@ -0,0 +1,130 @@
+#
+# 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__ = ["clean"]
+
+
+import numpy as np
+
+from lsst.meas.astrom.sip import LeastSqFitter1dPoly
+
+
+def clean(srcMatch, wcs, order=3, nsigma=3):
+    """Remove bad points from srcMatch
+
+    Input:
+    srcMatch : list of det::SourceMatch
+    order:      Order of polynomial to use in robust fitting
+    nsigma:    Sources more than this far away from the robust best fit
+                polynomial are removed
+
+    Return:
+    list of det::SourceMatch of the good data points
+    """
+
+    N = len(srcMatch)
+    catX = np.zeros(N)
+    # catY = np.zeros(N)
+    for i in range(N):
+        x, y = wcs.skyToPixel(srcMatch[i].first.getCoord())
+        catX[i] = x
+        # catY[i] = y
+
+    # TODO -- why does this only use X?
+
+    x = np.array([s.second.getX() for s in srcMatch])
+    dx = x - catX
+    sigma = np.zeros_like(dx) + 0.1
+
+    idx = indicesOfGoodPoints(x, dx, sigma, order=order, nsigma=nsigma)
+
+    clean = []
+    for i in idx:
+        clean.append(srcMatch[i])
+    return clean
+
+
+def indicesOfGoodPoints(x, y, s, order=1, nsigma=3, maxiter=100):
+    """Return a list of indices in the range [0, len(x)]
+    of points that lie less than nsigma away from the robust
+    best fit polynomial
+    """
+
+    # Indices of elements of x sorted in order of increasing value
+    idx = x.argsort()
+    newidx = []
+    for niter in range(maxiter):
+        rx = chooseRx(x, idx, order)
+        ry = chooseRy(y, idx, order)
+        rs = np.ones((len(rx)))
+
+        lsf = LeastSqFitter1dPoly(list(rx), list(ry), list(rs), order)
+        fit = [lsf.valueAt(value) for value in rx]
+        f = [lsf.valueAt(value) for value in x]
+
+        sigma = (y-f).std()
+        if sigma == 0:
+            # all remaining points are good; short circuit
+            return newidx if newidx else idx
+        deviance = np.fabs((y - f) / sigma)
+        newidx = np.flatnonzero(deviance < nsigma)
+
+        if False:
+            import matplotlib.pyplot as plt
+            plt.plot(x, y, 'ks')
+            plt.plot(rx, ry, 'b-')
+            plt.plot(rx, ry, 'bs')
+            plt.plot(rx, fit, 'ms')
+            plt.plot(rx, fit, 'm-')
+            # plt.plot(x[newidx], y[newidx], 'rs')
+            plt.show()
+
+        # If we haven't culled any points we're finished cleaning
+        if len(newidx) == len(idx):
+            break
+
+    # We get here because we either a) stopped finding bad points
+    # or b) ran out of iterations. Either way, we just return our
+    # list of indices of good points.
+    return newidx
+
+
+def chooseRx(x, idx, order):
+    """Create order+1 values of the ordinate based on the median of groups of elements of x"""
+    rSize = len(idx)/float(order+1)  # Note, a floating point number
+    rx = np.zeros((order+1))
+
+    for i in range(order+1):
+        rng = list(range(int(rSize*i), int(rSize*(i+1))))
+        rx[i] = np.mean(x[idx[rng]])
+    return rx
+
+
+def chooseRy(y, idx, order):
+    """Create order+1 values of the ordinate based on the median of groups of elements of y"""
+    rSize = len(idx)/float(order+1)  # Note, a floating point number
+    ry = np.zeros((order+1))
+
+    for i in range(order+1):
+        rng = list(range(int(rSize*i), int(rSize*(i+1))))
+        ry[i] = np.median(y[idx[rng]])
+    return ry

tests/test_anetAstrometryTask.py

@@ -0,0 +1,162 @@
+#
+# LSST Data Management System
+# Copyright 2008-2017 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/>.
+#
+
+import os.path
+import math
+import unittest
+
+import lsst.utils.tests
+import lsst.geom
+import lsst.afw.geom as afwGeom
+import lsst.afw.table as afwTable
+import lsst.afw.image as afwImage
+import lsst.meas.base as measBase
+from lsst.meas.extensions.astrometryNet import AstrometryNetDataConfig, \
+    ANetAstrometryTask, LoadAstrometryNetObjectsTask
+from test_findAstrometryNetDataDir import setupAstrometryNetDataDir
+
+
+class TestAstrometricSolver(lsst.utils.tests.TestCase):
+
+    def setUp(self):
+        self.datapath = setupAstrometryNetDataDir('photocal')
+
+        self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(3001, 3001))
+        crpix = lsst.geom.Box2D(self.bbox).getCenter()
+        self.tanWcs = afwGeom.makeSkyWcs(crpix=crpix,
+                                         crval=lsst.geom.SpherePoint(215.5, 53.0, lsst.geom.degrees),
+                                         cdMatrix=afwGeom.makeCdMatrix(scale=5.1e-5*lsst.geom.degrees))
+        self.exposure = afwImage.ExposureF(self.bbox)
+        self.exposure.setWcs(self.tanWcs)
+        self.exposure.setFilter(afwImage.Filter("r", True))
+        andConfig = AstrometryNetDataConfig()
+        andConfig.load(os.path.join(self.datapath, 'andConfig2.py'))
+        andConfig.magErrorColumnMap = {}
+        self.refObjLoader = LoadAstrometryNetObjectsTask(andConfig=andConfig)
+
+    def tearDown(self):
+        del self.tanWcs
+        del self.exposure
+        del self.refObjLoader
+
+    def testTrivial(self):
+        """Test fit with no distortion
+        """
+        self.doTest(afwGeom.makeIdentityTransform())
+
+    def testRadial(self):
+        """Test fit with radial distortion
+
+        The offset comes from the fact that the CCD is not centered
+        """
+        self.doTest(afwGeom.makeRadialTransform([0, 1.01, 1e-7]))
+
+    def makeSourceSchema(self):
+        schema = afwTable.SourceTable.makeMinimalSchema()
+        measBase.SingleFrameMeasurementTask(schema=schema)  # expand the schema
+        return schema
+
+    def doTest(self, pixelsToTanPixels):
+        """Test using pixelsToTanPixels to distort the source positions
+        """
+        distortedWcs = afwGeom.makeModifiedWcs(pixelTransform=pixelsToTanPixels, wcs=self.tanWcs,
+                                               modifyActualPixels=False)
+        self.exposure.setWcs(distortedWcs)
+        sourceCat = self.makeSourceCat(distortedWcs)
+        print("number of stars =", len(sourceCat))
+        config = ANetAstrometryTask.ConfigClass()
+        solver = ANetAstrometryTask(config=config, refObjLoader=self.refObjLoader,
+                                    schema=self.makeSourceSchema())
+        results = solver.run(
+            sourceCat=sourceCat,
+            exposure=self.exposure,
+        )
+        fitWcs = self.exposure.getWcs()
+        self.assertRaises(Exception, self.assertWcsAlmostEqualOverBBox, fitWcs, distortedWcs)
+        self.assertWcsAlmostEqualOverBBox(distortedWcs, fitWcs, self.bbox,
+                                          maxDiffSky=0.01*lsst.geom.arcseconds, maxDiffPix=0.02)
+
+        srcCoordKey = afwTable.CoordKey(sourceCat.schema["coord"])
+        refCoordKey = afwTable.CoordKey(results.refCat.schema["coord"])
+        maxAngSep = 0*lsst.geom.radians
+        maxPixSep = 0
+        for refObj, src, d in results.matches:
+            refCoord = refObj.get(refCoordKey)
+            refPixPos = fitWcs.skyToPixel(refCoord)
+            srcCoord = src.get(srcCoordKey)
+            srcPixPos = src.getCentroid()
+
+            angSep = refCoord.separation(srcCoord)
+            maxAngSep = max(maxAngSep, angSep)
+
+            pixSep = math.hypot(*(srcPixPos-refPixPos))
+            maxPixSep = max(maxPixSep, pixSep)
+        print("max angular separation = %0.4f arcsec" % (maxAngSep.asArcseconds(),))
+        print("max pixel separation = %0.3f" % (maxPixSep,))
+        self.assertLess(maxAngSep.asArcseconds(), 0.005)
+        self.assertLess(maxPixSep, 0.03)
+
+        # try again, but without fitting the WCS
+        config.forceKnownWcs = True
+        solverNoFit = ANetAstrometryTask(config=config, refObjLoader=self.refObjLoader,
+                                         schema=self.makeSourceSchema())
+        self.exposure.setWcs(distortedWcs)
+        noFitResults = solverNoFit.run(
+            sourceCat=sourceCat,
+            exposure=self.exposure,
+        )
+        self.assertGreater(len(noFitResults.refCat), 300)
+
+    def makeSourceCat(self, distortedWcs):
+        """Make a source catalog by reading the position reference stars and distorting the positions
+        """
+        loadRes = self.refObjLoader.loadPixelBox(bbox=self.bbox, wcs=distortedWcs, filterName="r")
+        refCat = loadRes.refCat
+        refCentroidKey = afwTable.Point2DKey(refCat.schema["centroid"])
+        refFluxRKey = refCat.schema["r_flux"].asKey()
+
+        sourceSchema = self.makeSourceSchema()
+        sourceCat = afwTable.SourceCatalog(sourceSchema)
+        sourceCentroidKey = afwTable.Point2DKey(sourceSchema["slot_Centroid"])
+        sourceFluxKey = sourceSchema["slot_PsfFlux_flux"].asKey()
+        sourceFluxSigmaKey = sourceSchema["slot_PsfFlux_fluxSigma"].asKey()
+
+        sourceCat.reserve(len(refCat))
+        for refObj in refCat:
+            src = sourceCat.addNew()
+            src.set(sourceCentroidKey, refObj.get(refCentroidKey))
+            src.set(sourceFluxKey, refObj.get(refFluxRKey))
+            src.set(sourceFluxSigmaKey, refObj.get(refFluxRKey)/100)
+        return sourceCat
+
+
+class MemoryTester(lsst.utils.tests.MemoryTestCase):
+    pass
+
+
+def setup_module(module):
+    lsst.utils.tests.init()
+
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()