DM-32055: Write util for calculating angle between wcses

python/lsst/afw/geom/skyWcs.py

@@ -20,8 +20,10 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
+import scipy
 
 from lsst.utils import continueClass
+import lsst.geom
 from ._python import reduceTransform
 from ._geom import (SkyWcs, makeCdMatrix, makeFlippedWcs, makeModifiedWcs,
                     makeSkyWcs, makeTanSipWcs, makeWcsPairTransform,
@@ -34,7 +36,7 @@
            "makeHpxWcs"]
 
 
-@continueClass  # noqa: F811 (FIXME: remove for py 3.8+)
+@continueClass
 class SkyWcs:  # noqa: F811
     def pixelToSkyArray(self, x, y, degrees=False):
         """
@@ -99,5 +101,50 @@ def skyToPixelArray(self, ra, dec, degrees=False):
 
         return x.ravel(), y.ravel()
 
+    def getRelativeRotationToWcs(self, otherWcs):
+        """Get the difference in sky rotation angle to the specified wcs.
+
+        Ignoring location on the sky, if another wcs were atop this one,
+        what would the difference in rotation be? i.e. for
+
+        otherWcs = createInitialSkyWcsFromBoresight(radec, rotation, detector)
+
+        what is the value that needs to be added to ``self.rotation`` (or
+        subtracted from `other.rotation``) to align them?
+
+        Parameters
+        ----------
+        otherWcs : `lsst.afw.geom.SkyWcs`
+            The wcs to calculate the angle to.
+
+        Returns
+        -------
+        angle : `lsst.geom.Angle`
+            The angle between this and the supplied wcs,
+            over the half-open range [0, 2pi).
+        """
+        # Note: tests for this function live in
+        # obs_lsst/tests/test_afwWcsUtil.py due to the need for an easy
+        # constructor and instantiated detector, and the fact that afw
+        # cannot depend on obs_base or obs_lsst.
+
+        m1 = self.getCdMatrix()
+        m2 = otherWcs.getCdMatrix()
+
+        svd1 = scipy.linalg.svd(m1)
+        svd2 = scipy.linalg.svd(m2)
+
+        m1rot = np.matmul(svd1[0], svd1[2])
+        m2rot = np.matmul(svd2[0], svd2[2])
+
+        v_rot = [1, 0]
+
+        v_rot = np.matmul(v_rot, m1rot)  # rotate by wcs1
+        v_rot = np.matmul(v_rot, m2rot.T)  # rotate _back_ by wcs2
+
+        rotation = np.arctan2(v_rot[1], v_rot[0])
+        rotation = rotation % (2*np.pi)
+        return lsst.geom.Angle(rotation, lsst.geom.radians)
+
 
 SkyWcs.__reduce__ = reduceTransform