Update for afw::image::Wcs->afw::geom::SkyWcs

include/lsst/meas/extensions/psfex/Field.hh

@@ -11,8 +11,8 @@ namespace lsst {
    namespace daf { namespace base {
       class PropertySet;
    }}
-   namespace afw { namespace image {
-      class Wcs;
+   namespace afw { namespace geom {
+      class SkyWcs;
    }}
 }
 
@@ -42,7 +42,7 @@ public:
     //
     void finalize() { _finalize(true); }        
 
-    void addExt(lsst::afw::image::Wcs const& wcs, int const naxis1, int const naxis2, int const nobj=0);
+    void addExt(lsst::afw::geom::SkyWcs const& wcs, int const naxis1, int const naxis2, int const nobj=0);
 
     /// Return the number of extensions
     int getNext() const { return impl->next; }

python/lsst/meas/extensions/psfex/field.cc

@@ -25,7 +25,7 @@
 
 #undef PI // defined by psfex, clashes with afw
 
-#include "lsst/afw/image/Wcs.h"
+#include "lsst/afw/geom/SkyWcs.h"
 
 namespace py = pybind11;
 using namespace pybind11::literals;

python/lsst/meas/extensions/psfex/psfex.py

@@ -16,8 +16,9 @@
 except ImportError:
     plt = None
 
-import lsst.afw.coord as afwCoord
+from lsst.afw.coord import IcrsCoord
 import lsst.afw.geom as afwGeom
+from lsst.afw.fits import readMetadata
 import lsst.afw.image as afwImage
 import lsst.afw.table as afwTable
 import lsst.afw.display.ds9 as ds9
@@ -596,7 +597,9 @@ def showPsf(psf, set, ext=None, wcsData=None, trim=0, nspot=5,
     for i in range(2):
         delta = pos[1][1][i].asDegrees() - pos[0][1][i].asDegrees()
         CD.append(delta/(pos[1][0][i] - pos[0][0][i]))
-    mosWcs = afwImage.makeWcs(pos[0][1], pos[0][0], CD[0], 0, 0, CD[1])
+    CD = np.array(CD)
+    CD.shape = (2, 2)
+    mosWcs = afwGeom.makeSkyWcs(crval=pos[0][0], crpix=pos[0][1], cdMatrix=CD)
 
     if ext is not None:
         title = "%s-%d" % (title, ext)
@@ -710,16 +713,19 @@ def makeitLsst(prefs, context, saveWcs=False, plot=dict()):
         with pyfits.open(cat):
             # Hack: I want the WCS so I'll guess where the calexp is to be found
             calexpFile = guessCalexp(cat)
-            md = afwImage.readMetadata(calexpFile)
-            wcs = afwImage.makeWcs(md)
+            md = readMetadata(calexpFile)
+            wcs = afwGeom.makeSkyWcs(md)
 
             if not wcs:
-                crval = afwCoord.makeCoord(afwCoord.ICRS, 0.0*afwGeom.degrees, 0.0*afwGeom.degrees)
-                wcs = afwImage.makeWcs(crval, afwGeom.PointD(0, 0), 1.0, 0, 0, 1.0)
+                cdMatrix = np.array([1.0, 0.0, 0.0, 1.0])
+                cdMatrix.shape = (2, 2)
+                wcs = afwGeom.makeSkyWcs(crpix=afwGeom.PointD(0, 0),
+                                         crval=IcrsCoord(0.0*afwGeom.degrees, 0.0*afwGeom.degrees),
+                                         cdMatrix=cdMatrix)
 
             naxis1, naxis2 = md.get("NAXIS1"), md.get("NAXIS2")
             # Find how many rows there are in the catalogue
-            md = afwImage.readMetadata(cat)
+            md = readMetadata(cat)
 
             field.addExt(wcs, naxis1, naxis2, md.get("NAXIS2"))
             if saveWcs:
@@ -1032,7 +1038,7 @@ def makeit(prefs, context, saveWcs=False, plot=dict()):
                         for k in md.names():
                             if re.search(r"A$", k):
                                 md.set(k[:-1], md.get(k))
-                    wcs = afwImage.makeWcs(md)
+                    wcs = afwGeom.makeSkyWcs(md)
                     naxis1, naxis2 = md.get("NAXIS1"), md.get("NAXIS2")
                 elif hdu.name == "LDAC_OBJECTS":
                     nobj = len(hdu.data)

src/fieldImpl.cc

@@ -1,9 +1,8 @@
 // -*- lsst-C++ -*-
 #include <cstring>
 #include "lsst/meas/extensions/psfex/Field.hh"
-#include "wcslib/wcs.h"
 #undef PI
-#include "lsst/afw/image/Wcs.h"
+#include "lsst/afw/geom/SkyWcs.h"
 
 extern "C" {
     #include "globals.h"
@@ -83,19 +82,8 @@ Field::getPsfs() const
     return _psfs;
 }
 
-/************************************************************************************************************/
-//
-// This class exists solely so that I can recover the protected data member _wcsInfo
-//
-namespace {
-struct PsfUnpack : private lsst::afw::image::Wcs {
-    PsfUnpack(lsst::afw::image::Wcs const& wcs) : Wcs(wcs) { }
-    const struct wcsprm* getWcsInfo() { return _wcsInfo; }
-};
-}
-
 void
-Field::addExt(lsst::afw::image::Wcs const& wcs_,
+Field::addExt(lsst::afw::geom::SkyWcs const& wcs_,
               int const naxis1, int const naxis2,
               int const nobj)
 {
@@ -106,44 +94,48 @@ Field::addExt(lsst::afw::image::Wcs const& wcs_,
     /*
      * We're going to fake psfex's wcsstruct object.  We only need enough of it for field_locate
      */
-    PsfUnpack wcsUnpacked(wcs_);
-    struct wcsprm const* wcsPrm = wcsUnpacked.getWcsInfo();
     QMALLOC(impl->wcs[impl->next], wcsstruct, 1);
     wcsstruct *wcs = impl->wcs[impl->next];
 
-    wcs->naxis = wcsPrm->naxis;
+    wcs->naxis = 2;
     wcs->naxisn[0] = naxis1;
     wcs->naxisn[1] = naxis2;
-
+
+    auto const crval = wcs_.getSkyOrigin();
+    // crpix using the FITS standard = pixel origin using the LSST standard + 1
+    auto const crpix = wcs_.getPixelOrigin() + afw::geom::Extent2D(1, 1);
+    auto const cdMatrix = wcs_.getCdMatrix();
+    std::string const cunit("DEG");
+    auto metadata = wcs_.getFitsMetadata(false);
     for (int i = 0; i != wcs->naxis; ++i) {
-        strncpy(wcs->ctype[i], wcsPrm->ctype[i], sizeof(wcs->ctype[i]) - 1);
-        strncpy(wcs->cunit[i], wcsPrm->cunit[i], sizeof(wcs->cunit[i]) - 1);
-        wcs->crval[i] = wcsPrm->crval[i];
-
-        wcs->cdelt[i] = wcsPrm->cdelt[i];
-        wcs->crpix[i] = wcsPrm->crpix[i];
-        wcs->crder[i] = wcsPrm->crder[i];
-        wcs->csyer[i] = wcsPrm->csyer[i];
-        wcs->crval[i] = wcsPrm->crval[i];
+        auto ifits = i + 1;
+        auto ctype = metadata->getAsString("CTYPE" + std::to_string(ifits));
+        strncpy(wcs->ctype[i], ctype.c_str(), ctype.size() + 1);
+        strncpy(wcs->cunit[i], cunit.c_str(), cunit.size() + 1);
+        wcs->crpix[i] = crpix[i];
+        wcs->crval[i] = crval[i].asDegrees();
+        wcs->cdelt[i] = 1.0;  // scale is in the CD matrix (is this even needed?)
+        wcs->crder[i] = 0;
+        wcs->csyer[i] = 0;
     }
-    for (int i = 0; i != wcs->naxis*wcs->naxis; ++i) {
-        wcs->cd[i] = wcsPrm->cd[i];
+    for (int i = 0, k = 0; i < 2; ++i) {
+        for (int j = 0; j < 2; ++j, ++k) {
+            wcs->cd[k] = cdMatrix(i, j);
+        }
     }
-    wcs->longpole = wcsPrm->lonpole;
-    wcs->latpole = wcsPrm->latpole;
-    wcs->lat = wcsPrm->lat;
-    wcs->lng = wcsPrm->lng;
-    wcs->equinox = wcsPrm->equinox;
+    wcs->lng = 0;
+    wcs->lat = 1;
+    wcs->equinox = 2000;
 
-    CONST_PTR(afw::coord::Coord) center = wcs_.pixelToSky(afw::geom::Point2D(0.5*naxis1, 0.5*naxis2));
-    wcs->wcsscalepos[0] = center->getLongitude().asDegrees();
-    wcs->wcsscalepos[1] = center->getLatitude().asDegrees();
+    auto center = wcs_.pixelToSky(afw::geom::Point2D(0.5*naxis1, 0.5*naxis2));
+    wcs->wcsscalepos[0] = center.getLongitude().asDegrees();
+    wcs->wcsscalepos[1] = center.getLatitude().asDegrees();
 
     double maxradius = 0.0;             // Maximum distance to wcsscalepos
     for (int x = 0; x <= 1; ++x) {
         for (int y = 0; y <= 1; ++y) {
             afw::geom::Point2D point(x*naxis1, y*naxis2); // Corner
-            double const radius = center->angularSeparation(*wcs_.pixelToSky(point)).asDegrees();
+            double const radius = center.angularSeparation(wcs_.pixelToSky(point)).asDegrees();
             if (radius > maxradius) {
                 maxradius = radius;
             }

tests/testPsfexPsf.py

@@ -9,7 +9,7 @@
 #
 # 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
+# the Free Software Foundation, either vers[1.0ion 3 of the License, or
 # (at your option) any later version.
 #
 # This program is distributed in the hope that it will be useful,
@@ -30,7 +30,7 @@
 
 import lsst.utils.tests
 import lsst.afw.image as afwImage
-import lsst.afw.coord as afwCoord
+from lsst.afw.coord import IcrsCoord
 import lsst.afw.detection as afwDetection
 import lsst.afw.geom as afwGeom
 import lsst.afw.math as afwMath
@@ -104,8 +104,11 @@ def setUp(self):
         self.exposure = afwImage.makeExposure(self.mi)
         self.exposure.setPsf(measAlg.DoubleGaussianPsf(self.ksize, self.ksize,
                                                        1.5*sigma1, 1, 0.1))
-        crval = afwCoord.makeCoord(afwCoord.ICRS, 0.0*afwGeom.degrees, 0.0*afwGeom.degrees)
-        wcs = afwImage.makeWcs(crval, afwGeom.PointD(0, 0), 1.0, 0, 0, 1.0)
+        cdMatrix = np.array([1.0, 0.0, 0.0, 1.0])
+        cdMatrix.shape = (2, 2)
+        wcs = afwGeom.makeSkyWcs(crpix=afwGeom.PointD(0, 0),
+                                 crval=IcrsCoord(0.0*afwGeom.degrees, 0.0*afwGeom.degrees),
+                                 cdMatrix=cdMatrix)
         self.exposure.setWcs(wcs)
 
         #