Merge pull request #703 from lsst/tickets/DM-15180

DM-15180: Add coordinate covariance to coordinate calculation

include/lsst/afw/table/Source.h

@@ -24,7 +24,6 @@
 #ifndef AFW_TABLE_Source_h_INCLUDED
 #define AFW_TABLE_Source_h_INCLUDED
 
-
 #include "lsst/afw/detection/Footprint.h"
 #include "lsst/afw/table/Simple.h"
 #include "lsst/afw/table/aggregates.h"
@@ -33,6 +32,7 @@
 #include "lsst/afw/table/BaseColumnView.h"
 #include "lsst/afw/table/slots.h"
 #include "lsst/afw/table/io/FitsWriter.h"
+#include "lsst/afw/table/wcsUtils.h"
 
 namespace lsst {
 namespace afw {
@@ -188,10 +188,10 @@ class SourceRecord : public SimpleRecord {
     double getIxy() const;
 
     /// Update the coord field using the given Wcs and the field in the centroid slot.
-    void updateCoord(geom::SkyWcs const &wcs);
+    void updateCoord(geom::SkyWcs const &wcs, bool include_covariance = true);
 
     /// Update the coord field using the given Wcs and the image center from the given key.
-    void updateCoord(geom::SkyWcs const &wcs, PointKey<double> const &key);
+    void updateCoord(geom::SkyWcs const &wcs, PointKey<double> const &key,  bool include_covariance = true);
 
     SourceRecord(const SourceRecord &) = delete;
     SourceRecord &operator=(const SourceRecord &) = delete;

include/lsst/afw/table/aggregates.h

@@ -42,6 +42,8 @@ class Quadrupole;
 
 namespace table {
 
+template <typename T, int N>
+class CovarianceMatrixKey;
 /**
  *  A FunctorKey used to get or set a lsst::geom::Point from an (x,y) pair of int or double Keys.
  */
@@ -298,6 +300,11 @@ class CoordKey : public FunctorKey<lsst::geom::SpherePoint> {
      *  @param[in]     doc     String used as the documentation for the fields.
      */
     static CoordKey addFields(afw::table::Schema& schema, std::string const& name, std::string const& doc);
+
+    using ErrorKey = CovarianceMatrixKey<float, 2>;
+    static ErrorKey getErrorKey(Schema const & schema);
+
+    static ErrorKey addErrorFields(Schema & schema);
 
     /// Default constructor; instance will not be usable unless subsequently assigned to.
     CoordKey() noexcept : _ra(), _dec() {}
@@ -624,6 +631,27 @@ class CovarianceMatrixKey : public FunctorKey<Eigen::Matrix<T, N, N> > {
     ErrKeyArray _err;
     CovarianceKeyArray _cov;
 };
+
+class CoordErrKey : public CovarianceMatrixKey<float, 2> {
+public:
+    using ErrKeyArray = std::vector<Key<float>>;
+    using CovarianceKeyArray = std::vector<Key<float>>;
+    using NameArray = std::vector<std::string>;
+
+    static CoordErrKey addFields(Schema& schema, std::string const& name, std::string const& unit,
+                                 bool diagonalOnly = false);
+    static CoordErrKey addFields(Schema& schema, std::string const& name, NameArray const& units,
+                                 bool diagonalOnly = false);
+
+    CoordErrKey();
+
+    explicit CoordErrKey(ErrKeyArray const& err,
+                         CovarianceKeyArray const& cov = CovarianceKeyArray());
+    CoordErrKey(SubSchema const& s, std::string const& name);
+private:
+    ErrKeyArray _err;
+    CovarianceKeyArray _cov;
+};
 }  // namespace table
 }  // namespace afw
 }  // namespace lsst

include/lsst/afw/table/wcsUtils.h

@@ -56,14 +56,27 @@ void updateRefCentroids(geom::SkyWcs const& wcs, ReferenceCollection& refList);
  * @param[in,out] sourceList  Collection of sources. The schema must have two fields:
  *                  - "slot_Centroid": a field containing lsst::geom::Point2D; this field is read
  *                  - "coord": a field containing an ICRS lsst::afw::SpherePoint; this field is written
- *
+ * @param[in] include_covariance  Update coordinate uncertainty using the centroid uncertainty.
+ * 
  * @throws lsst::pex::exceptions::NotFoundError if refList's schema does not have the required fields.
  */
 template <typename SourceCollection>
-void updateSourceCoords(geom::SkyWcs const& wcs, SourceCollection& sourceList);
+void updateSourceCoords(geom::SkyWcs const& wcs, SourceCollection& sourceList, bool include_covariance=true);
+
+/**
+ * Calculate covariance for sky coordinates
+ *
+ * 
+ * @param[in] wcs  WCS to map from pixels to sky.
+ * @param[in] center  The object centroid.
+ * @param[in] err  Covariance matrix of the centroid.
+ * 
+ */
+Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const& wcs, lsst::geom::Point2D center, Eigen::Matrix2f err);
 
 }  // namespace table
 }  // namespace afw
 }  // namespace lsst
 
+
 #endif  // #ifndef LSST_AFW_TABLE_WCSUTILS_H

python/lsst/afw/table/_aggregates.cc

@@ -138,6 +138,8 @@ static void declareCoordKey(WrapperCollection &wrappers) {
         cls.def("__eq__", &CoordKey::operator==, py::is_operator());
         cls.def("__ne__", &CoordKey::operator!=, py::is_operator());
         cls.def_static("addFields", &CoordKey::addFields, "schema"_a, "name"_a, "doc"_a);
+        cls.def_static("addErrorFields", &CoordKey::addErrorFields, "schema"_a);
+        cls.def_static("getErrorKey", &CoordKey::getErrorKey, "schema"_a);
         cls.def("getRa", &CoordKey::getRa);
         cls.def("getDec", &CoordKey::getDec);
         cls.def("isValid", &CoordKey::isValid);

python/lsst/afw/table/_source.cc

@@ -103,12 +103,12 @@ PySourceRecord declareSourceRecord(WrapperCollection &wrappers) {
         cls.def("getIxx", &SourceRecord::getIxx);
         cls.def("getIyy", &SourceRecord::getIyy);
         cls.def("getIxy", &SourceRecord::getIxy);
-        cls.def("updateCoord", (void (SourceRecord::*)(geom::SkyWcs const &)) & SourceRecord::updateCoord,
-                "wcs"_a);
+        cls.def("updateCoord", (void (SourceRecord::*)(geom::SkyWcs const &, bool)) & SourceRecord::updateCoord,
+                "wcs"_a, "include_covariance"_a=true);
         cls.def("updateCoord",
-                (void (SourceRecord::*)(geom::SkyWcs const &, PointKey<double> const &)) &
+                (void (SourceRecord::*)(geom::SkyWcs const &, PointKey<double> const &, bool)) &
                         SourceRecord::updateCoord,
-                "wcs"_a, "key"_a);
+                "wcs"_a, "key"_a, "include_covariance"_a=true);
     });
 }
 

python/lsst/afw/table/_wcsUtils.cc

@@ -55,7 +55,8 @@ void declareUpdateRefCentroids(WrapperCollection &wrappers) {
 template <typename SourceCollection>
 void declareUpdateSourceCoords(WrapperCollection &wrappers) {
     wrappers.wrap([](auto &mod) {
-        mod.def("updateSourceCoords", updateSourceCoords<SourceCollection>, "wcs"_a, "sourceList"_a);
+        mod.def("updateSourceCoords", updateSourceCoords<SourceCollection>, "wcs"_a, "sourceList"_a,
+                "include_covariance"_a=true);
     });
 }
 

src/table/Source.cc

@@ -365,10 +365,28 @@ static SourceFitsReader const sourceFitsReader;
 
 SourceRecord::~SourceRecord() = default;
 
-void SourceRecord::updateCoord(geom::SkyWcs const &wcs) { setCoord(wcs.pixelToSky(getCentroid())); }
+void SourceRecord::updateCoord(geom::SkyWcs const &wcs, bool include_covariance) {
+    lsst::geom::Point2D center = getCentroid();
+    setCoord(wcs.pixelToSky(center));
+    if (include_covariance) {
+        // Get coordinate covariance:
+        auto err = getCentroidErr();
+        CoordKey::ErrorKey const coordErrKey = CoordKey::getErrorKey(getSchema());
+        Eigen::Matrix2f skyCov = calculateCoordCovariance(wcs, center, err);
+        set(coordErrKey, skyCov);
+    }
+}
 
-void SourceRecord::updateCoord(geom::SkyWcs const &wcs, PointKey<double> const &key) {
-    setCoord(wcs.pixelToSky(get(key)));
+void SourceRecord::updateCoord(geom::SkyWcs const &wcs, PointKey<double> const &key, bool include_covariance) {
+    lsst::geom::Point2D center = get(key);
+    setCoord(wcs.pixelToSky(center));
+    if (include_covariance) {
+        // Get coordinate covariance:
+        auto err = getCentroidErr();
+        CoordKey::ErrorKey const coordErrKey = CoordKey::getErrorKey(getSchema());
+        Eigen::Matrix2f skyCov = calculateCoordCovariance(wcs, center, err);
+        set(coordErrKey, skyCov);
+    }
 }
 
 void SourceRecord::_assign(BaseRecord const &other) {

src/table/aggregates.cc

@@ -121,6 +121,14 @@ void CoordKey::set(BaseRecord &record, lsst::geom::SpherePoint const &value) con
     record.set(_dec, value.getLatitude());
 }
 
+CoordKey::ErrorKey CoordKey::getErrorKey(Schema const & schema) {
+    return ErrorKey(schema["coord"], {"ra", "dec"});
+}
+
+CoordKey::ErrorKey CoordKey::addErrorFields(Schema & schema) {
+    return ErrorKey::addFields(schema, "coord", {"ra", "dec"}, "rad");
+}
+
 //============ QuadrupoleKey ================================================================================
 
 QuadrupoleKey QuadrupoleKey::addFields(Schema &schema, std::string const &name, std::string const &doc,

src/table/wcsUtils.cc

@@ -91,33 +91,81 @@ void updateRefCentroids(geom::SkyWcs const &wcs, ReferenceCollection &refList) {
     }
 }
 
+Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const& wcs, lsst::geom::Point2D center, Eigen::Matrix2f err) {
+    // Get the derivative of the pixel-to-sky transformation, then use it to
+    // propagate the centroid uncertainty to coordinate uncertainty. Note that
+    // the calculation is done in arcseconds, then converted to radians in
+    // order to achieve higher precision.
+    double scale = 1.0 / 3600.0;
+    Eigen::Matrix2d cdMatrix { {scale, 0}, {0, scale}};
+
+    lsst::geom::SpherePoint skyCenter = wcs.pixelToSky(center);
+
+    auto localGnomonicWcs = geom::makeSkyWcs(center, skyCenter, cdMatrix);
+    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::Matrix2f skyCov = d * err * d.transpose();
+
+    // Multiply by declination correction matrix in order to get sigma(RA) * cos(Dec) for the uncertainty
+    // in RA, and cov(RA, Dec) * cos(Dec) for the RA/Dec covariance:
+    float cosDec = std::cos(skyCenter.getDec().asRadians());
+    Eigen::Matrix2f decCorr{ {cosDec, 0}, {0, 1.0}};
+    Eigen::Matrix2f skyCovCorr = decCorr * skyCov * decCorr;
+    return skyCovCorr;
+}
+
+
 template <typename SourceCollection>
-void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList) {
+void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList, bool include_covariance) {
     if (sourceList.empty()) {
         return;
     }
     auto const schema = getSchema(sourceList[0]);
     Point2DKey const centroidKey(schema["slot_Centroid"]);
+    CovarianceMatrixKey<float, 2> const centroidErrKey(schema["slot_Centroid"], {"x", "y"});
     CoordKey const coordKey(schema["coord"]);
     std::vector<lsst::geom::Point2D> pixelList;
     pixelList.reserve(sourceList.size());
+    std::vector<Eigen::Matrix2f> skyErrList;
+    skyErrList.reserve(sourceList.size());
+
     for (auto const &source : sourceList) {
-        pixelList.emplace_back(getValue(source, centroidKey));
+        lsst::geom::Point2D center = getValue(source, centroidKey);
+        pixelList.emplace_back(center);
+        if (include_covariance) {
+            auto err = getValue(source, centroidErrKey);
+            Eigen::Matrix2f skyCov = calculateCoordCovariance(wcs, center, err);
+            skyErrList.emplace_back(skyCov);
+        }
     }
+
     std::vector<lsst::geom::SpherePoint> const skyList = wcs.pixelToSky(pixelList);
     auto skyCoord = skyList.cbegin();
+    auto skyErr = skyErrList.cbegin();
     for (auto &source : sourceList) {
         setValue(source, coordKey, *skyCoord);
+        if (include_covariance) {
+            CoordKey::ErrorKey const coordErrKey = CoordKey::getErrorKey(schema);
+            setValue(source, coordErrKey, *skyErr);
+        }
         ++skyCoord;
+        ++skyErr;
     }
 }
 
 /// @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>> &);
-template void updateSourceCoords(geom::SkyWcs const &, SourceCatalog &);
+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);
 /// @endcond
 
 }  // namespace table

tests/test_sourceTable.py

@@ -112,6 +112,8 @@ def setUp(self):
             self.schema, "d", "", lsst.afw.table.CoordinateType.PIXEL)
         self.psfShapeFlagKey = self.schema.addField("d_flag", type="Flag")
 
+        self.coordErrKey = lsst.afw.table.CoordKey.addErrorFields(self.schema)
+
         self.table = lsst.afw.table.SourceTable.make(self.schema)
         self.catalog = lsst.afw.table.SourceCatalog(self.table)
         self.record = self.catalog.addNew()
@@ -223,6 +225,31 @@ def testCoordUpdate(self):
         coord2 = wcs.pixelToSky(self.record.get(self.centroidKey))
         self.assertEqual(coord1, coord2)
 
+    def testCoordErrors(self):
+        self.table.defineCentroid("b")
+        wcs = makeWcs()
+        self.record.updateCoord(wcs)
+
+        scale = (1.0 * lsst.geom.arcseconds).asDegrees()
+        center = self.record.getCentroid()
+        skyCenter = wcs.pixelToSky(center)
+        localGnomonicWcs = lsst.afw.geom.makeSkyWcs(
+            center, skyCenter, np.diag((scale, scale)))
+        measurementToLocalGnomonic = wcs.getTransform().then(
+            localGnomonicWcs.getTransform().inverted()
+        )
+        localMatrix = measurementToLocalGnomonic.getJacobian(center)
+        radMatrix = np.radians(localMatrix / 3600)
+
+        centroidErr = self.record.getCentroidErr()
+        coordErr = radMatrix.dot(centroidErr.dot(radMatrix.T))
+        # multiply RA by cos(Dec):
+        cosDec = np.cos(skyCenter.getDec().asRadians())
+        decCorr = np.array([[cosDec, 0], [0, 1]])
+        coordErrCorr = decCorr.dot(coordErr.dot(decCorr))
+        catCoordErr = self.record.get(self.coordErrKey)
+        np.testing.assert_almost_equal(coordErrCorr, catCoordErr, decimal=16)
+
     def testSorting(self):
         self.assertFalse(self.catalog.isSorted())
         self.catalog.sort()

tests/test_tableUtils.py

@@ -64,6 +64,10 @@ def setUp(self):
         srcSchema = afwTable.SourceTable.makeMinimalSchema()
         self.srcCentroidKey = afwTable.Point2DKey.addFields(srcSchema, "base_SdssCentroid",
                                                             "centroid", "pixels")
+        self.srcCentroidErrKey = afwTable.CovarianceMatrix2fKey.addFields(srcSchema, "base_SdssCentroid",
+                                                                          ["x", "y"], "pixels")
+        self.srcCoordErrKey = afwTable.CoordKey.addErrorFields(srcSchema)
+
         srcAliases = srcSchema.getAliasMap()
         srcAliases.set("slot_Centroid", "base_SdssCentroid")
         self.srcCoordKey = afwTable.CoordKey(srcSchema["coord"])
@@ -146,6 +150,33 @@ def testCatalogs(self):
         # check that centroids and coords match
         self.checkCatalogs()
 
+    def testCoordErrors(self):
+        """Check that updateSourceCoords has correctly propagated the centroid
+        errors.
+        """
+        maxPix = 2000
+        numPoints = 10
+
+        self.setCatalogs(maxPix=maxPix, numPoints=numPoints)
+        scale = (1.0 * lsst.geom.arcseconds).asDegrees()
+        # update the catalogs
+        afwTable.updateSourceCoords(self.wcs, self.sourceCat)
+        for src in self.sourceCat:
+            center = src.get(self.srcCentroidKey)
+            skyCenter = self.wcs.pixelToSky(center)
+            localGnomonicWcs = lsst.afw.geom.makeSkyWcs(
+                center, skyCenter, np.diag((scale, scale)))
+            measurementToLocalGnomonic = self.wcs.getTransform().then(
+                localGnomonicWcs.getTransform().inverted()
+            )
+            localMatrix = measurementToLocalGnomonic.getJacobian(center)
+            radMatrix = np.radians(localMatrix / 3600)
+
+            centroidErr = src.get(self.srcCentroidErrKey)
+            coordErr = radMatrix.dot(centroidErr.dot(radMatrix.T))
+            catCoordErr = src.get(self.srcCoordErrKey)
+            np.testing.assert_almost_equal(coordErr, catCoordErr)
+
     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"""
         self.assertEqual(len(self.sourceCat), len(self.refCat))
@@ -182,8 +213,10 @@ def setCatalogs(self, maxPix, numPoints):
         for i in np.linspace(-maxPix, maxPix, numPoints):
             for j in np.linspace(-maxPix, maxPix, numPoints):
                 centroid = lsst.geom.Point2D(i, j)
+                centroidErr = np.array([[0.05, 0], [0, 0.05]])
                 src = self.sourceCat.addNew()
                 src.set(self.srcCentroidKey, centroid)
+                src.set(self.srcCentroidErrKey, centroidErr)
 
                 refObj = self.refCat.addNew()
                 coord = self.wcs.pixelToSky(centroid)

tests/test_ticketDM-433.py

@@ -135,6 +135,7 @@ def setUp(self):
         self.makeInstFlux(self.schema, "a", 1)
         self.makeCentroid(self.schema, "b", 2)
         self.makeShape(self.schema, "c", 2)
+        lsst.afw.table.CoordKey.addErrorFields(self.schema)
         self.table = lsst.afw.table.SourceTable.make(self.schema)
         self.catalog = lsst.afw.table.SourceCatalog(self.table)
         self.record = self.catalog.addNew()