From 1b6092986d6e7d969fc3a71d5bbcb9512254577c Mon Sep 17 00:00:00 2001 From: Clare Saunders Date: Mon, 19 Jun 2023 12:14:47 -0700 Subject: [PATCH] Add coordinate covariance and keys to access it --- include/lsst/afw/table/Source.h | 6 +-- include/lsst/afw/table/aggregates.h | 28 ++++++++++++++ include/lsst/afw/table/wcsUtils.h | 17 ++++++++- python/lsst/afw/table/_aggregates.cc | 2 + python/lsst/afw/table/_source.cc | 8 ++-- python/lsst/afw/table/_wcsUtils.cc | 3 +- src/table/Source.cc | 24 ++++++++++-- src/table/aggregates.cc | 8 ++++ src/table/wcsUtils.cc | 56 ++++++++++++++++++++++++++-- tests/test_sourceTable.py | 27 ++++++++++++++ tests/test_tableUtils.py | 33 ++++++++++++++++ tests/test_ticketDM-433.py | 1 + 12 files changed, 196 insertions(+), 17 deletions(-) diff --git a/include/lsst/afw/table/Source.h b/include/lsst/afw/table/Source.h index 412d244bf6..47d0e01db8 100644 --- a/include/lsst/afw/table/Source.h +++ b/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 const &key); + void updateCoord(geom::SkyWcs const &wcs, PointKey const &key, bool include_covariance = true); SourceRecord(const SourceRecord &) = delete; SourceRecord &operator=(const SourceRecord &) = delete; diff --git a/include/lsst/afw/table/aggregates.h b/include/lsst/afw/table/aggregates.h index 34ead15b5a..bc0939ae98 100644 --- a/include/lsst/afw/table/aggregates.h +++ b/include/lsst/afw/table/aggregates.h @@ -42,6 +42,8 @@ class Quadrupole; namespace table { +template +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 { * @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; + 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 > { ErrKeyArray _err; CovarianceKeyArray _cov; }; + +class CoordErrKey : public CovarianceMatrixKey { +public: + using ErrKeyArray = std::vector>; + using CovarianceKeyArray = std::vector>; + using NameArray = std::vector; + + 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 diff --git a/include/lsst/afw/table/wcsUtils.h b/include/lsst/afw/table/wcsUtils.h index 6d6a8eeafb..47d2a7b5b6 100644 --- a/include/lsst/afw/table/wcsUtils.h +++ b/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 -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 diff --git a/python/lsst/afw/table/_aggregates.cc b/python/lsst/afw/table/_aggregates.cc index 9fc8101f2c..e5b023ed85 100644 --- a/python/lsst/afw/table/_aggregates.cc +++ b/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); diff --git a/python/lsst/afw/table/_source.cc b/python/lsst/afw/table/_source.cc index 53a0ddbdd8..8506e54898 100644 --- a/python/lsst/afw/table/_source.cc +++ b/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 const &)) & + (void (SourceRecord::*)(geom::SkyWcs const &, PointKey const &, bool)) & SourceRecord::updateCoord, - "wcs"_a, "key"_a); + "wcs"_a, "key"_a, "include_covariance"_a=true); }); } diff --git a/python/lsst/afw/table/_wcsUtils.cc b/python/lsst/afw/table/_wcsUtils.cc index 31392bed90..2f8a90e876 100644 --- a/python/lsst/afw/table/_wcsUtils.cc +++ b/python/lsst/afw/table/_wcsUtils.cc @@ -55,7 +55,8 @@ void declareUpdateRefCentroids(WrapperCollection &wrappers) { template void declareUpdateSourceCoords(WrapperCollection &wrappers) { wrappers.wrap([](auto &mod) { - mod.def("updateSourceCoords", updateSourceCoords, "wcs"_a, "sourceList"_a); + mod.def("updateSourceCoords", updateSourceCoords, "wcs"_a, "sourceList"_a, + "include_covariance"_a=true); }); } diff --git a/src/table/Source.cc b/src/table/Source.cc index 8b9d7a1efa..82bef4b35c 100644 --- a/src/table/Source.cc +++ b/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 const &key) { - setCoord(wcs.pixelToSky(get(key))); +void SourceRecord::updateCoord(geom::SkyWcs const &wcs, PointKey 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) { diff --git a/src/table/aggregates.cc b/src/table/aggregates.cc index 7e78e0b311..b0864054e5 100644 --- a/src/table/aggregates.cc +++ b/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, diff --git a/src/table/wcsUtils.cc b/src/table/wcsUtils.cc index 375bdcdbf4..4426a93a0e 100644 --- a/src/table/wcsUtils.cc +++ b/src/table/wcsUtils.cc @@ -91,24 +91,70 @@ 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 -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 const centroidErrKey(schema["slot_Centroid"], {"x", "y"}); CoordKey const coordKey(schema["coord"]); std::vector pixelList; pixelList.reserve(sourceList.size()); + std::vector 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 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; } } @@ -116,8 +162,10 @@ void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList) { template void updateRefCentroids(geom::SkyWcs const &, std::vector> &); template void updateRefCentroids(geom::SkyWcs const &, SimpleCatalog &); -template void updateSourceCoords(geom::SkyWcs const &, std::vector> &); -template void updateSourceCoords(geom::SkyWcs const &, SourceCatalog &); +template void updateSourceCoords(geom::SkyWcs const &, std::vector> &, + bool include_covariance); +template void updateSourceCoords(geom::SkyWcs const &, SourceCatalog &, + bool include_covariance); /// @endcond } // namespace table diff --git a/tests/test_sourceTable.py b/tests/test_sourceTable.py index 7cbfc3d2ec..57aae693ec 100644 --- a/tests/test_sourceTable.py +++ b/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() diff --git a/tests/test_tableUtils.py b/tests/test_tableUtils.py index 87d1e303eb..94c2a78d25 100644 --- a/tests/test_tableUtils.py +++ b/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) diff --git a/tests/test_ticketDM-433.py b/tests/test_ticketDM-433.py index 5b068ad942..c4d98ef767 100644 --- a/tests/test_ticketDM-433.py +++ b/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()