DM-10973: Make SkyWcs transform to IcrsCoord instead of SpherePoint

include/lsst/afw/geom/Endpoint.h

@@ -30,8 +30,8 @@
 #include "astshim.h"
 #include "ndarray.h"
 
+#include "lsst/afw/coord/Coord.h"
 #include "lsst/afw/geom/Point.h"
-#include "lsst/afw/geom/SpherePoint.h"
 
 namespace lsst {
 namespace afw {
@@ -41,10 +41,10 @@ namespace geom {
 Virtual base class for endpoints, which are helper classes for Transform
 
 Endpoints transform points and lists of points from LSST-specific data types,
-such as Point2D and SpherePoint, to a form accepted by ast::Mapping.tran.
+such as Point2D and IcrsCoord, to a form accepted by ast::Mapping.tran.
 Each type of endpoint is used for a particular LSST data type, for example:
 - Point2Endpoint is used for Point2D data
-- SpherePointEndpoint for SpherePoint data
+- IcrsCoordEndpoint for IcrsCoord data
 - GenericEndpoint is used when no other form will do; its LSST data type
   is identical to the type used for ast::Mapping.applyForward.
 
@@ -164,7 +164,7 @@ class BaseEndpoint {
 /**
 Base class for endpoints with Array = std::vector<Point> where Point has 2 dimensions
 
-@note Subclasses must provide `arrayFromData` 
+@note Subclasses must provide `dataFromPoint`, `dataFromArray`, `pointFromData` and `arrayFromData`
 */
 template <typename PointT>
 class BaseVectorEndpoint : public BaseEndpoint<PointT, std::vector<PointT>> {
@@ -179,13 +179,7 @@ class BaseVectorEndpoint : public BaseEndpoint<PointT, std::vector<PointT>> {
 
     virtual ~BaseVectorEndpoint(){};
 
-    virtual int getNPoints(Array const &arr) const override { return arr.size(); }
-
-    virtual std::vector<double> dataFromPoint(Point const &point) const override;
-
-    virtual ndarray::Array<double, 2, 2> dataFromArray(Array const &arr) const override;
-
-    virtual Point pointFromData(std::vector<double> const &data) const override;
+    virtual int getNPoints(Array const &arr) const override;
 
 protected:
     /**
@@ -262,6 +256,12 @@ class Point2Endpoint : public BaseVectorEndpoint<Point2D> {
 
     virtual ~Point2Endpoint(){};
 
+    virtual std::vector<double> dataFromPoint(Point const &point) const override;
+
+    virtual ndarray::Array<double, 2, 2> dataFromArray(Array const &arr) const override;
+
+    virtual Point pointFromData(std::vector<double> const &data) const override;
+
     virtual Array arrayFromData(ndarray::Array<double, 2, 2> const &data) const override;
 
     /**
@@ -277,24 +277,24 @@ class Point2Endpoint : public BaseVectorEndpoint<Point2D> {
 };
 
 /**
-An endpoint for SpherePoint
+An endpoint for IcrsCoord
 
-A SpherePointEndpoint always has 2 axes: longitude, latitude
+A IcrsCoordEndpoint always has 2 axes: RA, Dec
 */
-class SpherePointEndpoint : public BaseVectorEndpoint<SpherePoint> {
+class IcrsCoordEndpoint : public BaseVectorEndpoint<coord::IcrsCoord> {
 public:
-    SpherePointEndpoint(SpherePointEndpoint const &) = default;
-    SpherePointEndpoint(SpherePointEndpoint &&) = default;
-    SpherePointEndpoint &operator=(SpherePointEndpoint const &) = delete;
-    SpherePointEndpoint &operator=(SpherePointEndpoint &&) = delete;
+    IcrsCoordEndpoint(IcrsCoordEndpoint const &) = default;
+    IcrsCoordEndpoint(IcrsCoordEndpoint &&) = default;
+    IcrsCoordEndpoint &operator=(IcrsCoordEndpoint const &) = delete;
+    IcrsCoordEndpoint &operator=(IcrsCoordEndpoint &&) = delete;
 
     /**
-    Construct a SpherePointEndpoint
+    Construct a IcrsCoordEndpoint
     */
-    explicit SpherePointEndpoint() : BaseVectorEndpoint(2) {}
+    explicit IcrsCoordEndpoint() : BaseVectorEndpoint(2) {}
 
     /**
-    Construct a SpherePointEndpoint with nAxes specified; nAxes must equal 2
+    Construct a IcrsCoordEndpoint with nAxes specified; nAxes must equal 2
 
     This constructor is primarily used by Transform; other users are encouraged
     to use the default constructor.
@@ -303,9 +303,15 @@ class SpherePointEndpoint : public BaseVectorEndpoint<SpherePoint> {
 
     @throws lsst.pex.exceptions.InvalidParameterError if nAxes != 2
     */
-    explicit SpherePointEndpoint(int nAxes);
+    explicit IcrsCoordEndpoint(int nAxes);
 
-    virtual ~SpherePointEndpoint(){};
+    virtual ~IcrsCoordEndpoint(){};
+
+    virtual std::vector<double> dataFromPoint(Point const &point) const override;
+
+    virtual ndarray::Array<double, 2, 2> dataFromArray(Array const &arr) const override;
+
+    virtual Point pointFromData(std::vector<double> const &data) const override;
 
     virtual Array arrayFromData(ndarray::Array<double, 2, 2> const &data) const override;
 
@@ -328,11 +334,11 @@ std::ostream &operator<<(std::ostream &os, GenericEndpoint const &endpoint);
 /// Print "Point2Endpoint()" to the ostream
 std::ostream &operator<<(std::ostream &os, Point2Endpoint const &endpoint);
 
-/// Print "SpherePointEndpoint()" to the ostream
-std::ostream &operator<<(std::ostream &os, SpherePointEndpoint const &endpoint);
+/// Print "IcrsCoordEndpoint()" to the ostream
+std::ostream &operator<<(std::ostream &os, IcrsCoordEndpoint const &endpoint);
 
-}  // geom
-}  // afw
-}  // lsst
+}  // namespace geom
+}  // namespace afw
+}  // namespace lsst
 
 #endif

include/lsst/afw/geom/SkyWcs.h

@@ -32,9 +32,9 @@
 #include "astshim.h"
 #include "ndarray.h"
 
+#include "lsst/afw/coord/Coord.h"
 #include "lsst/afw/geom/Angle.h"
 #include "lsst/afw/geom/Point.h"
-#include "lsst/afw/geom/SpherePoint.h"
 #include "lsst/afw/geom/Endpoint.h"
 #include "lsst/afw/geom/Transform.h"
 #include "lsst/daf/base/PropertyList.h"
@@ -94,7 +94,7 @@ Constructors must set the following properties of the contained ast::FrameSet:
     SkyWcs uses this to look up CRVAL; "Ignored" is required to prevent CRVAL from being used as an offset
     from the desired WCS.
 */
-class SkyWcs : public Transform<Point2Endpoint, SpherePointEndpoint> {
+class SkyWcs : public Transform<Point2Endpoint, IcrsCoordEndpoint> {
 public:
     SkyWcs(SkyWcs const &) = default;
     SkyWcs(SkyWcs &&) = default;
@@ -110,7 +110,7 @@ class SkyWcs : public Transform<Point2Endpoint, SpherePointEndpoint> {
     @param[in] cdMatrix  CD matrix, where element (i-1, j-1) corresponds to FITS keyword CDi_j
                         and i, j have range [1, 2]. May be computed by calling makeCdMatrix.
     */
-    explicit SkyWcs(Point2D const &crpix, SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix);
+    explicit SkyWcs(Point2D const &crpix, coord::IcrsCoord const &crval, Eigen::Matrix2d const &cdMatrix);
 
     /**
     Construct a WCS from FITS keywords
@@ -160,7 +160,7 @@ class SkyWcs : public Transform<Point2Endpoint, SpherePointEndpoint> {
     /**
     Get CRVAL, the sky origin or celestial fiducial point
     */
-    SpherePoint getSkyOrigin() const;
+    coord::IcrsCoord getSkyOrigin() const;
 
     /**
     Get the 2x2 CD matrix at the specified pixel position
@@ -189,8 +189,8 @@ class SkyWcs : public Transform<Point2Endpoint, SpherePointEndpoint> {
     */
     //@{
     std::pair<Angle, Angle> pixelToSky(double x, double y) const;
-    SpherePoint pixelToSky(Point2D const &pixel) const { return applyForward(pixel); };
-    std::vector<SpherePoint> pixelToSky(std::vector<Point2D> const &pixels) const {
+    coord::IcrsCoord pixelToSky(Point2D const &pixel) const { return applyForward(pixel); };
+    std::vector<coord::IcrsCoord> pixelToSky(std::vector<Point2D> const &pixels) const {
         return applyForward(pixels);
     }
     //@}
@@ -202,8 +202,10 @@ class SkyWcs : public Transform<Point2Endpoint, SpherePointEndpoint> {
     */
     //@{
     std::pair<double, double> skyToPixel(Angle const &ra, Angle const &dec) const;
-    Point2D skyToPixel(SpherePoint const &sky) const { return applyInverse(sky); }
-    std::vector<Point2D> skyToPixel(std::vector<SpherePoint> const &sky) const { return applyInverse(sky); }
+    Point2D skyToPixel(coord::IcrsCoord const &sky) const { return applyInverse(sky); }
+    std::vector<Point2D> skyToPixel(std::vector<coord::IcrsCoord> const &sky) const {
+        return applyInverse(sky);
+    }
     //@}
 
 private:

include/lsst/afw/geom/Transform.h

@@ -38,7 +38,7 @@ namespace afw {
 namespace geom {
 
 /**
-Transform LSST spatial data, such as Point2D and SpherePoint, using an AST transform.
+Transform LSST spatial data, such as Point2D and IcrsCoord, using an AST transform.
 
 This class contains two Endpoints, to specify the "from" and "to" LSST data type,
 and an ast::FrameSet or ast::Mapping to specify the transformation.
@@ -216,8 +216,7 @@ class Transform {
      *     auto pixelsToSky = pixelsToFP.then(fpToPupil).then(pupilToSky);
      */
     template <class NextToEndpoint>
-    Transform<FromEndpoint, NextToEndpoint> then(
-            Transform<ToEndpoint, NextToEndpoint> const &next) const;
+    Transform<FromEndpoint, NextToEndpoint> then(Transform<ToEndpoint, NextToEndpoint> const &next) const;
 
 protected:
     /**
@@ -243,8 +242,8 @@ for example "Transform<GenericEndpoint(4), Point2Endpoint()>"
 template <class FromEndpoint, class ToEndpoint>
 std::ostream &operator<<(std::ostream &os, Transform<FromEndpoint, ToEndpoint> const &transform);
 
-}  // geom
-}  // afw
-}  // lsst
+}  // namespace geom
+}  // namespace afw
+}  // namespace lsst
 
 #endif

include/lsst/afw/geom/detail/frameSetUtils.h

@@ -30,6 +30,7 @@
 #include "astshim.h"
 #include "ndarray.h"
 
+#include "lsst/afw/coord/Coord.h"
 #include "lsst/afw/geom/Endpoint.h"
 #include "lsst/afw/geom/Transform.h"
 #include "lsst/daf/base/PropertyList.h"
@@ -60,7 +61,8 @@ Make FITS metadata for a pure tangent WCS
 @param[in] cdMatrix  CD matrix where element (i-1, j-1) corresponds to FITS keyword CDi_j
                         and i, j have range [1, 2]
 */
-std::shared_ptr<daf::base::PropertyList> makeTanWcsMetadata(Point2D const& crpix, SpherePoint const& crval,
+std::shared_ptr<daf::base::PropertyList> makeTanWcsMetadata(Point2D const& crpix,
+                                                            coord::IcrsCoord const& crval,
                                                             Eigen::Matrix2d const& cdMatrix);
 
 /**

include/lsst/afw/geom/transformFactory.h

@@ -55,23 +55,23 @@ Transform<FromEndpoint, ToEndpoint> linearizeTransform(
         typename Transform<FromEndpoint, ToEndpoint>::FromPoint const &point);
 
 /*
- * The correct behavior for linearization is unclear where SpherePoints are involved (see discussion on
+ * The correct behavior for linearization is unclear where IcrsCoords are involved (see discussion on
  * DM-10542). Forbid usage until somebody needs it. Note to maintainers: the template specializations MUST
  * be deleted in the header for compilers to complain correctly.
  */
 #define DISABLE(From, To)                                                         \
     template <>                                                                   \
     Transform<From, To> linearizeTransform<From, To>(Transform<From, To> const &, \
                                                      Transform<From, To>::FromPoint const &) = delete;
-DISABLE(GenericEndpoint, SpherePointEndpoint);
-DISABLE(Point2Endpoint, SpherePointEndpoint);
-DISABLE(SpherePointEndpoint, GenericEndpoint);
-DISABLE(SpherePointEndpoint, Point2Endpoint);
-DISABLE(SpherePointEndpoint, SpherePointEndpoint);
+DISABLE(GenericEndpoint, IcrsCoordEndpoint);
+DISABLE(Point2Endpoint, IcrsCoordEndpoint);
+DISABLE(IcrsCoordEndpoint, GenericEndpoint);
+DISABLE(IcrsCoordEndpoint, Point2Endpoint);
+DISABLE(IcrsCoordEndpoint, IcrsCoordEndpoint);
 #undef DISABLE
 
-}  // geom
-}  // afw
-}  // lsst
+}  // namespace geom
+}  // namespace afw
+}  // namespace lsst
 
 #endif  // LSST_AFW_GEOM_TRANSFORMFACTORY_H

python/lsst/afw/coord/utils.py

@@ -28,6 +28,8 @@
 
 import warnings
 
+import numpy as np
+
 import lsst.utils.tests
 import lsst.afw.geom as afwGeom
 
@@ -54,6 +56,31 @@ def assertCoordsAlmostEqual(testCase, coord0, coord1, maxDiff=0.001*afwGeom.arcs
                       (msg, measDiff.asArcseconds(), maxDiff.asArcseconds()))
 
 
+@lsst.utils.tests.inTestCase
+def assertCoordListsAlmostEqual(testCase, coordlist0, coordlist1, maxDiff=0.001*afwGeom.arcseconds, msg=None):
+    """!Assert that two lists of IcrsCoords represent almost the same point on the sky
+
+    @warning the coordinate systems are not compared; instead both sets of angles are converted to ICRS
+    and the angular separation measured.
+
+    @param[in] testCase  unittest.TestCase instance the test is part of;
+                        an object supporting one method: fail(self, msgStr)
+    @param[in] coordlist0  list of IcrsCoords 0
+    @param[in] coordlist1  list of IcrsCoords 1
+    @param[in] maxDiff  maximum angular separation, an lsst.afw.geom.Angle
+    @param[in] msg  exception message prefix; details of the error are appended after ": "
+    """
+    testCase.assertEqual(len(coordlist0), len(coordlist1), msg=msg)
+    sepArr = np.array([sp0.toIcrs().angularSeparation(sp1.toIcrs())
+                       for sp0, sp1 in zip(coordlist0, coordlist1)])
+    badArr = sepArr > maxDiff
+    if np.any(badArr):
+        maxInd = np.argmax(sepArr)
+        testCase.fail("%s: IcrsCoordLists differ in %s places; max separation is at %s: %s\" > %s\"%s" %
+                      (msg, np.sum(badArr), maxInd, sepArr[maxInd].asArcseconds(),
+                       maxDiff.asArcseconds()))
+
+
 @lsst.utils.tests.inTestCase
 def assertCoordsNearlyEqual(*args, **kwargs):
     warnings.warn("Deprecated. Use assertCoordsAlmostEqual",

python/lsst/afw/geom/detail/frameSetUtils.cc

@@ -30,7 +30,6 @@
 #include "lsst/daf/base.h"
 #include "lsst/afw/geom/Angle.h"
 #include "lsst/afw/geom/Point.h"
-#include "lsst/afw/geom/SpherePoint.h"
 #include "lsst/afw/geom/detail/frameSetUtils.h"
 
 namespace py = pybind11;