DM-10765: Replace existing WCS classes with SkyWcs

examples/ticket2710.py

@@ -9,7 +9,8 @@
 
 import lsst.meas.astrom as measAstrom
 import lsst.afw.table as afwTable
-import lsst.afw.image as afwImage
+from lsst.afw.fits import readMetadata
+from lsst.afw.geom import makeSkyWcs
 from lsst.log import Log
 
 import pyfits
@@ -157,10 +158,8 @@ def showSipSolutions(srcs, wcs0, andDir, x0, y0, W, H, filterName,
         msy.append(src.getY())
 
     plt.clf()
-    #plt.plot(x, y, 'o', mec='r', mfc='none', ms=4)
     plt.plot(x, y, 'r.')
     plt.plot(msx, msy, 'o', mec='r')
-    #plt.plot(rx0, ry0, 'o', mec='g', mfc='none', ms=4)
     plt.plot(rx0, ry0, 'g.')
     plt.plot(mrx, mry, 'gx')
     plt.title('TAN matches')
@@ -170,7 +169,6 @@ def showSipSolutions(srcs, wcs0, andDir, x0, y0, W, H, filterName,
 
     solve = ast.determineWcs2(srcs, **imargs)
     wcs1 = solve.sipWcs
-    # print 'wcs1:', wcs1.getFitsMetadata().toString()
 
     matches = solve.sipMatches
     msx, msy = [], []
@@ -249,6 +247,7 @@ def readSourcesFromXyTable(xyfn):
         src.set(fkey, fi)
     return srcs
 
+
 if __name__ == '__main__':
     mydir = os.path.dirname(__file__)
     # fitscopy coaddSources.fits"[col x=centroid_sdss[1]; y=centroid_sdss[2]; flux=flux_psf]" xy.fits
@@ -261,8 +260,8 @@ def readSourcesFromXyTable(xyfn):
 
     # Read original WCS
     fn = os.path.join(mydir, 't2710.wcs')
-    hdr = afwImage.readMetadata(fn)
-    wcs0 = afwImage.makeWcs(hdr)
+    hdr = readMetadata(fn)
+    wcs0 = makeSkyWcs(hdr)
 
     anddir = os.path.join(mydir, 'astrometry_net_data', 'ticket2710')
 

include/lsst/meas/extensions/astrometryNet/astrometry_net.h

@@ -35,8 +35,6 @@ extern "C" {
 #include "lsst/daf/base/Persistable.h"
 #include "lsst/daf/base/PropertyList.h"
 #include "lsst/afw/table.h"
-#include "lsst/afw/image/Wcs.h"
-#include "lsst/afw/image/TanWcs.h"
 #include "lsst/afw/geom.h"
 #include "lsst/afw/coord.h"
 
@@ -167,7 +165,7 @@ class Solver {
 
     std::shared_ptr<lsst::daf::base::PropertyList> getSolveStats() const;
 
-    std::shared_ptr<lsst::afw::image::Wcs> getWcs();
+    std::shared_ptr<lsst::afw::geom::SkyWcs> getWcs();
 
     bool didSolve() const {
         return solver_did_solve(_solver.get());

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

@@ -26,7 +26,6 @@
 from builtins import input
 from builtins import zip
 from builtins import range
-from contextlib import contextmanager
 
 import numpy as np
 
@@ -287,42 +286,6 @@ def distort(self, sourceCat, exposure):
 
         return afwGeom.Box2I(bboxD)
 
-    @contextmanager
-    def distortionContext(self, sourceCat, exposure):
-        """!Context manager that applies and removes distortion
-
-        We move the "centroid" definition in the catalog table to
-        point to the distorted positions.  This is undone on exit
-        from the context.
-
-        The input Wcs is taken to refer to the coordinate system
-        with the distortion correction applied, and hence no shift
-        is required when the sources are distorted.  However, after
-        Wcs fitting, the Wcs is in the distorted frame so when the
-        distortion correction is removed, the Wcs needs to be
-        shifted to compensate.
-
-        \param sourceCat Sources on exposure, an lsst.afw.table.SourceCatalog
-        \param exposure Exposure holding Wcs, an lsst.afw.image.ExposureF or D
-        \return bounding box of distorted exposure
-        """
-        # Apply distortion, if not already present in the exposure's WCS
-        if exposure.getWcs().hasDistortion():
-            yield exposure.getBBox()
-        else:
-            bbox = self.distort(sourceCat=sourceCat, exposure=exposure)
-            oldCentroidName = sourceCat.table.getCentroidDefinition()
-            sourceCat.table.defineCentroid(self.distortedName)
-            try:
-                yield bbox  # Execute 'with' block, providing bbox to 'as' variable
-            finally:
-                # Un-apply distortion
-                sourceCat.table.defineCentroid(oldCentroidName)
-                x0, y0 = exposure.getXY0()
-                wcs = exposure.getWcs()
-                if wcs:
-                    wcs.shiftReferencePixel(-bbox.getMinX() + x0, -bbox.getMinY() + y0)
-
     @pipeBase.timeMethod
     def loadAndMatch(self, exposure, sourceCat, bbox=None):
         """!Load reference objects overlapping an exposure and match to sources detected on that exposure
@@ -339,36 +302,36 @@ def loadAndMatch(self, exposure, sourceCat, bbox=None):
 
         @note ignores config.forceKnownWcs
         """
-        with self.distortionContext(sourceCat=sourceCat, exposure=exposure) as bbox:
-            if not self.solver:
-                self.makeSubtask("solver")
-
-            astrom = self.solver.useKnownWcs(
-                sourceCat=sourceCat,
-                exposure=exposure,
-                bbox=bbox,
-                calculateSip=False,
-            )
-
-            if astrom is None or astrom.getWcs() is None:
-                raise RuntimeError("Unable to solve astrometry")
-
-            matches = astrom.getMatches()
-            matchMeta = astrom.getMatchMetadata()
-            if matches is None or len(matches) == 0:
-                raise RuntimeError("No astrometric matches")
-            self.log.info("%d astrometric matches" % (len(matches)))
-
-            if self._display:
-                frame = lsstDebug.Info(__name__).frame
-                displayAstrometry(exposure=exposure, sourceCat=sourceCat, matches=matches,
-                                  frame=frame, pause=False)
-
-            return pipeBase.Struct(
-                refCat=astrom.refCat,
-                matches=matches,
-                matchMeta=matchMeta,
-            )
+        bbox = exposure.getBBox()
+        if not self.solver:
+            self.makeSubtask("solver")
+
+        astrom = self.solver.useKnownWcs(
+            sourceCat=sourceCat,
+            exposure=exposure,
+            bbox=bbox,
+            calculateSip=False,
+        )
+
+        if astrom is None or astrom.getWcs() is None:
+            raise RuntimeError("Unable to solve astrometry")
+
+        matches = astrom.getMatches()
+        matchMeta = astrom.getMatchMetadata()
+        if matches is None or len(matches) == 0:
+            raise RuntimeError("No astrometric matches")
+        self.log.info("%d astrometric matches" % (len(matches)))
+
+        if self._display:
+            frame = lsstDebug.Info(__name__).frame
+            displayAstrometry(exposure=exposure, sourceCat=sourceCat, matches=matches,
+                              frame=frame, pause=False)
+
+        return pipeBase.Struct(
+            refCat=astrom.refCat,
+            matches=matches,
+            matchMeta=matchMeta,
+        )
 
     @pipeBase.timeMethod
     def _astrometry(self, sourceCat, exposure, bbox=None):

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

@@ -529,7 +529,7 @@ def getBlindWcsSolution(self, sourceCat,
         if wcs is not None:
             if pixelScale is None:
                 if usePixelScale:
-                    pixelScale = wcs.pixelScale()
+                    pixelScale = wcs.getPixelScale()
                     self.log.debug('Setting pixel scale estimate = %.3f from given WCS estimate',
                                    pixelScale.asArcseconds())
 
@@ -550,7 +550,7 @@ def getBlindWcsSolution(self, sourceCat,
                                    'image size, and searchRadiusScale = %g',
                                    searchRadius, searchRadiusScale)
             if useParity:
-                parity = wcs.isFlipped()
+                parity = wcs.isFlipped
                 self.log.debug('Using parity = %s' % (parity and 'True' or 'False'))
 
         if doTrim:
@@ -612,7 +612,7 @@ def getSipWcsFromWcs(self, wcs, bbox, ngrid=20, linearizeAtCenter=True):
                 src.set(key.getX(), x0 + xx)
                 src.set(key.getY(), y0 + yy)
                 csrc.append(src)
-                rd = wcs.pixelToSky(afwGeom.Point2D(xx + x0, yy + y0))
+                rd = wcs.pixelToSky(xx + x0, yy + y0)
                 ref = refs.makeRecord()
                 ref.setCoord(rd)
                 cref.append(ref)
@@ -681,6 +681,11 @@ def _calculateSipTerms(self, origWcs, refCat, sourceCat, matches, bbox):
                 self.log.warn('Failed to calculate distortion terms. Error: ', str(e))
                 break
 
+            # update the source catalog
+            for source in sourceCat:
+                skyPos = proposedWcs.pixelToSky(source.getCentroid())
+                source.setCoord(skyPos)
+
             # use new WCS to get new matchlist.
             proposedMatchlist = self._getMatchList(sourceCat, refCat, proposedWcs)
             proposedMatchSize = len(proposedMatchlist)
@@ -950,12 +955,9 @@ def _solve(self, sourceCat, wcs, bbox, pixelScale, radecCenter, searchRadius, pa
         if solver.didSolve():
             self.log.debug('Solved!')
             wcs = solver.getWcs()
-            self.log.debug('WCS: %s', wcs.getFitsMetadata().toString())
 
             if x0 != 0 or y0 != 0:
-                wcs.shiftReferencePixel(x0, y0)
-                self.log.debug('After shifting reference pixel by x0,y0 = (%i,%i), WCS is: %s',
-                               x0, y0, wcs.getFitsMetadata().toString())
+                wcs = wcs.copyAtShiftedPixelOrigin(afwGeom.Extent2D(x0, y0))
 
         else:
             self.log.warn('Did not get an astrometric solution from Astrometry.net')
@@ -1009,11 +1011,11 @@ def _createMetadata(bbox, wcs, filterName):
 
     bboxD = afwGeom.Box2D(bbox)
     cx, cy = bboxD.getCenter()
-    radec = wcs.pixelToSky(cx, cy).toIcrs()
+    radec = wcs.pixelToSky(cx, cy)
     meta.add('RA', radec.getRa().asDegrees(), 'field center in degrees')
     meta.add('DEC', radec.getDec().asDegrees(), 'field center in degrees')
     pixelRadius = math.hypot(*bboxD.getDimensions())/2.0
-    skyRadius = wcs.pixelScale() * pixelRadius
+    skyRadius = wcs.getPixelScale() * pixelRadius
     meta.add('RADIUS', skyRadius.asDegrees(),
              'field radius in degrees, approximate')
     meta.add('SMATCHV', 1, 'SourceMatchVector version number')

src/astrometry_net.cc

@@ -138,19 +138,22 @@ std::shared_ptr<lsst::daf::base::PropertyList> Solver::getSolveStats() const {
     return qa;
 }
 
-std::shared_ptr<lsst::afw::image::Wcs> Solver::getWcs() {
+std::shared_ptr<lsst::afw::geom::SkyWcs> Solver::getWcs() {
     MatchObj* match = solver_get_best_match(_solver.get());
     if (!match)
-        return std::shared_ptr<afw::image::Wcs>();
+        return std::shared_ptr<afw::geom::SkyWcs>();
     tan_t* wcs = &(match->wcstan);
 
     afw::geom::Point2D crpix(wcs->crpix[0], wcs->crpix[1]);
-    std::shared_ptr<afw::coord::Coord const> crval
-        (new afw::coord::Coord(wcs->crval[0] * afw::geom::degrees,
-                             wcs->crval[1] * afw::geom::degrees));
-    return afw::image::makeWcs(*crval, crpix,
-                             wcs->cd[0][0], wcs->cd[0][1],
-                             wcs->cd[1][0], wcs->cd[1][1]);
+    auto const crval = afw::coord::IcrsCoord(wcs->crval[0] * afw::geom::degrees,
+                                             wcs->crval[1] * afw::geom::degrees);
+    Eigen::Matrix2d cdMatrix;
+    for (int i = 0; i < 2; ++i) {
+        for (int j = 0; j < 2; ++j) {
+            cdMatrix(i, j) = wcs->cd[i][j];
+        }
+    }
+    return afw::geom::makeSkyWcs(crpix, crval, cdMatrix);
 }
 
 void Solver::run(double cpulimit) {

tests/test_createWcsWithSip.py

@@ -74,7 +74,6 @@ def testBigXy0(self):
         res = self.astrom.determineWcs2(cat, bbox=bbox)
         self.assertIsNotNone(res.sipWcs, "Failed to fit SIP terms")
         print('Got result', res)
-        print('SIP:', res.sipWcs.getFitsMetadata().toString())
 
         wcs = res.wcs
         for src in cat:
@@ -106,30 +105,14 @@ def singleTestInstance(self, filename, distortFunc):
         res = self.astrom.determineWcs2(img, bbox=bbox)
         imgWcs = res.getWcs()
 
-        def printWcs(wcs):
-            print("WCS metadata:")
-            md = wcs.getFitsMetadata()
-            for name in md.names():
-                print("%s: %r" % (name, md.get(name)))
-
-        printWcs(imgWcs)
-
         # Create a wcs with sip
         cat = cat.cast(SimpleCatalog, False)
         matchList = self.matchSrcAndCatalogue(cat, img, imgWcs)
         print("*** num matches =", len(matchList))
         return
         sipObject = sip.makeCreateWcsWithSip(matchList, imgWcs, 3)
 
-        # print 'Put in TAN Wcs:'
-        # print imgWcs.getFitsMetadata().toString()
         imgWcs = sipObject.getNewWcs()
-        # print 'Got SIP Wcs:'
-        # print imgWcs.getFitsMetadata().toString()
-
-        # Write out the SIP header
-        # afwImage.fits_write_imageF('createWcsWithSip.sip', afwImage.ImageF(0,0),
-        # imgWcs.getFitsMetadata())
 
         print('number of matches:', len(matchList), sipObject.getNPoints())
         scatter = sipObject.getScatterOnSky().asArcseconds()

tests/test_loadAstrometryNetObjects.py

@@ -29,7 +29,6 @@
 import lsst.utils.tests
 from lsst.daf.base import PropertySet
 import lsst.afw.geom as afwGeom
-from lsst.afw.image import makeWcs
 from lsst.afw.table import CoordKey, Point2DKey
 from lsst.meas.extensions.astrometryNet import LoadAstrometryNetObjectsTask, \
     AstrometryNetDataConfig
@@ -63,7 +62,7 @@ def setUp(self):
         metadata.set("CD1_2", 0.0)
         metadata.set("CD2_2", -5.1e-05)
         metadata.set("CD2_1", 0.0)
-        self.wcs = makeWcs(metadata)
+        self.wcs = afwGeom.makeSkyWcs(metadata)
         self.desNumStarsInPixelBox = 270
         self.desNumStarsInSkyCircle = 410