Merge pull request #687 from lsst/tickets/DM-39053

DM-39053: Add WCS metadata stripping code that doesn't create a WCS.

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

@@ -97,6 +97,13 @@ All frames are instances of ast::Frame except the SKY frame. All have 2 axes.
 */
 std::shared_ptr<ast::FrameDict> readLsstSkyWcs(daf::base::PropertySet& metadata, bool strip = true);
 
+/**
+Strip all WCS metadata from a header.
+
+@param[in,out] metadata  FITS header cards
+*/
+void stripWcsMetadata(daf::base::PropertySet &metadata);
+
 /**
 Copy values from an AST FitsChan into a PropertyList
 

include/lsst/afw/geom/wcsUtils.h

@@ -184,6 +184,8 @@ std::shared_ptr<daf::base::PropertyList> makeTanSipMetadata(
         Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA, Eigen::MatrixXd const& sipB,
         Eigen::MatrixXd const& sipAp, Eigen::MatrixXd const& sipBp);
 
+void stripWcsMetadata(daf::base::PropertySet &metadata);
+
 }  // namespace geom
 }  // namespace afw
 }  // namespace lsst

python/lsst/afw/geom/_wcsUtils.cc

@@ -67,6 +67,7 @@ void wrapWcsUtils(lsst::utils::python::WrapperCollection &wrappers) {
                         Eigen::MatrixXd const &, Eigen::MatrixXd const &, Eigen::MatrixXd const &,
                         Eigen::MatrixXd const &))makeTanSipMetadata,
                 "crpix"_a, "crval"_a, "cdMatrix"_a, "sipA"_a, "sipB"_a, "sipAp"_a, "sipBp"_a);
+        mod.def("stripWcsMetadata", stripWcsMetadata, "metadata"_a);
     });
 }
 

src/geom/detail/frameSetUtils.cc

@@ -131,6 +131,51 @@ std::shared_ptr<ast::FrameSet> readFitsWcs(daf::base::PropertySet& metadata, boo
     return frameSet;
 }
 
+void stripWcsMetadata(daf::base::PropertySet &metadata) {
+    // Exclude WCS A keywords because LSST uses them to store XY0
+    auto wcsANames = createTrivialWcsMetadata("A", lsst::geom::Point2I(0, 0))->names();
+    std::set<std::string> excludeNames(wcsANames.begin(), wcsANames.end());
+    // Ignore NAXIS1, NAXIS2 because if they are zero then AST will fail to read a WCS
+    // Ignore LTV1/2 because LSST adds it and this code should ignore it and not strip it
+    // Exclude comments and history to reduce clutter
+    std::set<std::string> moreNames{"NAXIS1", "NAXIS2", "LTV1", "LTV2", "COMMENT", "HISTORY"};
+    excludeNames.insert(moreNames.begin(), moreNames.end());
+
+    // Replace RADECSYS with RADESYS if only the former is present
+    if (metadata.exists("RADECSYS") && !metadata.exists("RADESYS")) {
+        metadata.set("RADESYS", metadata.getAsString("RADECSYS"));
+        metadata.remove("RADECSYS");
+    }
+
+    auto channel = getFitsChanFromPropertyList(metadata, excludeNames,
+                                               "Encoding=FITS-WCS, IWC=1, SipReplace=0, ReportLevel=3");
+    auto const initialNames = setFromVector(channel.getAllCardNames());
+
+    std::shared_ptr<ast::Object> obj;
+    try {
+        obj = channel.read();
+    } catch (std::runtime_error const&) {
+        return;
+    }
+
+    auto const finalNames = setFromVector(channel.getAllCardNames());
+
+    // FITS keywords that FitsChan stripped
+    std::set<std::string> namesChannelStripped;
+    std::set_difference(initialNames.begin(), initialNames.end(), finalNames.begin(), finalNames.end(),
+                        std::inserter(namesChannelStripped, namesChannelStripped.begin()));
+
+    // FITS keywords that FitsChan may strip that we want to keep in `metadata`
+    std::set<std::string> const namesToKeep = {"DATE-OBS", "MJD-OBS"};
+
+    std::set<std::string> namesToStrip;  // names to strip from metadata
+    std::set_difference(namesChannelStripped.begin(), namesChannelStripped.end(), namesToKeep.begin(),
+                        namesToKeep.end(), std::inserter(namesToStrip, namesToStrip.begin()));
+    for (auto const& name : namesToStrip) {
+        metadata.remove(name);
+    }
+}
+
 std::shared_ptr<ast::FrameDict> readLsstSkyWcs(daf::base::PropertySet& metadata, bool strip) {
     // Record CRPIX in GRID coordinates
     // so we can compute CRVAL after standardizing the SkyFrame to ICRS

src/geom/wcsUtils.cc

@@ -212,6 +212,10 @@ std::shared_ptr<daf::base::PropertyList> makeTanSipMetadata(
     return metadata;
 }
 
+void stripWcsMetadata(daf::base::PropertySet &metadata) {
+    detail::stripWcsMetadata(metadata);
+}
+
 }  // namespace geom
 }  // namespace afw
 }  // namespace lsst

tests/test_skyWcs.py

@@ -18,7 +18,8 @@
     TransformPoint2ToPoint2, TransformPoint2ToSpherePoint, makeRadialTransform, \
     SkyWcs, makeSkyWcs, makeCdMatrix, makeWcsPairTransform, \
     makeFlippedWcs, makeModifiedWcs, makeTanSipWcs, \
-    getIntermediateWorldCoordsToSky, getPixelToIntermediateWorldCoords
+    getIntermediateWorldCoordsToSky, getPixelToIntermediateWorldCoords, \
+    stripWcsMetadata
 from lsst.afw.geom import getCdMatrixFromMetadata, getSipMatrixFromMetadata, makeSimpleWcsMetadata
 from lsst.afw.geom.testUtils import makeSipIwcToPixel, makeSipPixelToIwc
 from lsst.afw.fits import makeLimitedFitsHeader
@@ -694,6 +695,15 @@ def testBasics(self):
         makeSkyWcs(self.metadata, strip=True)
         self.assertEqual(len(self.metadata.names(False)), 0)
 
+    def testBasicsStrip(self):
+        stripWcsMetadata(self.metadata)
+        self.assertEqual(len(self.metadata.names(False)), 0)
+        # The metadata should be unchanged if we attempt to strip it again
+        metadataCopy = self.metadata.deepCopy()
+        stripWcsMetadata(self.metadata)
+        for key in self.metadata.keys():
+            self.assertEqual(self.metadata[key], metadataCopy[key])
+
     def testNormalizationFk5(self):
         """Test that readLsstSkyWcs correctly normalizes FK5 1975 to ICRS
         """