DM-53940: Add convertCentroid for X/Y to RA/dec covariance

.gitignore

@@ -181,6 +181,7 @@ tests/data/kernel*.boost
 ups/*.cfgc
 extensions/
 .vscode/
+.idea/
 
 # Output of tests
 doc/htmlDir

include/lsst/afw/table/wcsUtils.h

@@ -64,15 +64,46 @@ template <typename SourceCollection>
 void updateSourceCoords(geom::SkyWcs const& wcs, SourceCollection& sourceList, bool include_covariance=true);
 
 /**
- * Calculate covariance for sky coordinates
+ * Calculate covariance for sky coordinates given a pixel centroid and errors.
+ *
+ * @tparam t Floating point type of the input errors/output covariance.
  *
- * 
  * @param[in] wcs  WCS to map from pixels to sky.
- * @param[in] center  The object centroid.
+ * @param[in] center  The object centroid in pixels.
  * @param[in] err  Covariance matrix of the centroid.
- * 
+ * @param[in] factor Factor to multiply the WCS matrix terms, which are in degrees.
+ *                   Default is pi/180 to convert to radians.
+ *
+ * @return The RA/Dec sky covariance matrix in degrees, multiplied by factor.
+ */
+template <typename t>
+Eigen::Matrix<t, 2, 2> calculateCoordCovariance(
+    geom::SkyWcs const& wcs, lsst::geom::Point2D center, Eigen::Matrix<t, 2, 2> err,
+    double factor=lsst::geom::PI/180.0
+);
+
+/**
+ * Convert an x/y centroid with errors into RA/dec.
+ *
+ * @param[in] wcs  WCS to map from pixels to sky.
+ * @param[in] x  The x centroid in pixels.
+ * @param[in] y  The y centroid in pixels.
+ * @param[in] xErr  The standard error on the x centroid in pixels.
+ * @param[in] yErr  The standard error on the y centroid in pixels.
+ * @param[in] xy_covariance  The xy covariance in pixels squared.
+ *
+ * @return A pair of tuples with the RA/dec centroid in degrees,
+ *         and the RA error, dec error and RA/dec covariance in degrees
+ *         (squared for the covariance).
  */
-Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const& wcs, lsst::geom::Point2D center, Eigen::Matrix2f err);
+std::pair<std::tuple<double, double>, std::tuple<double, double, double>> convertCentroid(
+    geom::SkyWcs const &wcs,
+    double x,
+    double y,
+    double xErr,
+    double yErr,
+    double xy_covariance=0
+);
 
 }  // namespace table
 }  // namespace afw

python/lsst/afw/table/_wcsUtils.cc

@@ -60,9 +60,18 @@ void declareUpdateSourceCoords(WrapperCollection &wrappers) {
     });
 }
 
+void declareConvertCentroid(WrapperCollection &wrappers) {
+    wrappers.wrap([](auto &mod) {
+        mod.def("convertCentroid", convertCentroid, "wcs"_a, "x"_a, "y"_a, "xErr"_a, "yErr"_a,
+                "xy_covariance"_a=0.);
+    });
+}
+
 }  // namespace
 
 void wrapWcsUtils(WrapperCollection &wrappers) {
+    declareConvertCentroid(wrappers);
+
     declareUpdateRefCentroids<std::vector<std::shared_ptr<lsst::afw::table::SimpleRecord>>>(wrappers);
     declareUpdateRefCentroids<lsst::afw::table::SimpleCatalog>(wrappers);
 

src/table/wcsUtils.cc

@@ -92,10 +92,11 @@ void updateRefCentroids(geom::SkyWcs const &wcs, ReferenceCollection &refList) {
     }
 }
 
-Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const &wcs, lsst::geom::Point2D center,
-                                         Eigen::Matrix2f err) {
+template <typename t>
+Eigen::Matrix<t, 2, 2> calculateCoordCovariance(geom::SkyWcs const &wcs, lsst::geom::Point2D center,
+                                                Eigen::Matrix<t, 2, 2> err, double factor) {
     if (!std::isfinite(center.getX()) || !std::isfinite(center.getY())) {
-        return Eigen::Matrix2f::Constant(NAN);
+        return Eigen::Matrix<t, 2, 2>::Constant(NAN);
     }
     // Get the derivative of the pixel-to-sky transformation, then use it to
     // propagate the centroid uncertainty to coordinate uncertainty. Note that
@@ -109,9 +110,9 @@ Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const &wcs, lsst::geom::Po
     auto measurementToLocalGnomonic = wcs.getTransform()->then(*localGnomonicWcs->getTransform()->inverted());
 
     Eigen::Matrix2d localMatrix = measurementToLocalGnomonic->getJacobian(center);
-    Eigen::Matrix2f d = localMatrix.cast<float>() * scale * (lsst::geom::PI / 180.0);
+    Eigen::Matrix<t, 2, 2>  d = localMatrix.cast<t>() * scale * factor;
 
-    Eigen::Matrix2f skyCov = d * err * d.transpose();
+    Eigen::Matrix<t, 2, 2> skyCov = d * err * d.transpose();
     // Because the transform from pixels to RA/Dec was done at a local
     // gnomonic, the RA and Dec covariance are already commensurate, and
     // multiplying the RA covariance by cos(Dec) is not necessary.
@@ -137,7 +138,7 @@ void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList, b
         pixelList.emplace_back(center);
         if (include_covariance) {
             auto err = getValue(source, centroidErrKey);
-            Eigen::Matrix2f skyCov = calculateCoordCovariance(wcs, center, err);
+            Eigen::Matrix2f skyCov = calculateCoordCovariance<float>(wcs, center, err);
             skyErrList.emplace_back(skyCov);
         }
     }
@@ -156,13 +157,42 @@ void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList, b
     }
 }
 
+std::pair<std::tuple<double, double>, std::tuple<double, double, double>> convertCentroid(
+    geom::SkyWcs const &wcs,
+    double x,
+    double y,
+    double xErr,
+    double yErr,
+    double xy_covariance
+) {
+    lsst::geom::Point2D center{x, y};
+
+    auto radec = wcs.pixelToSky(center);
+
+    Eigen::Matrix2d err;
+    err(0, 0) = xErr*xErr;
+    err(1, 1) = yErr*yErr;
+    err(0, 1) = xy_covariance;
+    err(1, 0) = xy_covariance;
+
+    Eigen::Matrix2d result = calculateCoordCovariance<double>(wcs, center, err, 1.0);
+
+    return {
+        {radec.getRa().asDegrees(), radec.getDec().asDegrees()},
+        {sqrt(result(0, 0)), sqrt(result(1, 1)), result(0, 1)}
+    };
+}
+
 /// @cond
 template void updateRefCentroids(geom::SkyWcs const &, std::vector<std::shared_ptr<SimpleRecord>> &);
 template void updateRefCentroids(geom::SkyWcs const &, SimpleCatalog &);
 
 template void updateSourceCoords(geom::SkyWcs const &, std::vector<std::shared_ptr<SourceRecord>> &,
                                  bool include_covariance);
 template void updateSourceCoords(geom::SkyWcs const &, SourceCatalog &, bool include_covariance);
+
+template Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const &wcs, lsst::geom::Point2D center, Eigen::Matrix2f err, double factor);
+template Eigen::Matrix2d calculateCoordCovariance(geom::SkyWcs const &wcs, lsst::geom::Point2D center, Eigen::Matrix2d err, double factor);
 /// @endcond
 
 }  // namespace table

tests/test_tableUtils.py

@@ -161,6 +161,8 @@ def testCoordErrors(self):
         scale = (1.0 * lsst.geom.arcseconds).asDegrees()
         # update the catalogs
         afwTable.updateSourceCoords(self.wcs, self.sourceCat)
+        factor = math.pi/180.
+        factor_sq = factor**2
         for src in self.sourceCat:
             center = src.get(self.srcCentroidKey)
             skyCenter = self.wcs.pixelToSky(center)
@@ -175,7 +177,23 @@ def testCoordErrors(self):
             centroidErr = src.get(self.srcCentroidErrKey)
             coordErr = radMatrix.dot(centroidErr.dot(radMatrix.T))
             catCoordErr = src.get(self.srcCoordErrKey)
-            np.testing.assert_almost_equal(coordErr, catCoordErr)
+            np.testing.assert_allclose(coordErr, catCoordErr, rtol=1e-10, atol=1e-10)
+
+            (ra, dec), (ra_err, dec_err, ra_dec_cov) = afwTable.convertCentroid(
+                self.wcs,
+                center.getX(), center.getY(),
+                math.sqrt(centroidErr[0, 0]), math.sqrt(centroidErr[1, 1]), centroidErr[0, 1],
+            )
+            np.testing.assert_allclose(
+                (ra, dec),
+                (skyCenter[0].asDegrees(), skyCenter[1].asDegrees()),
+                rtol=1e-10, atol=1e-10,
+            )
+            np.testing.assert_allclose(
+                ((ra_err*factor)**2, (dec_err*factor)**2, ra_dec_cov*factor_sq),
+                (catCoordErr[0, 0], catCoordErr[1, 1], catCoordErr[0, 1]),
+                rtol=1e-10, atol=1e-10,
+            )
 
     def checkCatalogs(self, maxPixDiff=1e-5, maxSkyDiff=0.001*lsst.geom.arcseconds):
         """Check that the source and reference object catalogs have equal centroids and coords"""