Merge branch 'develop'

doc/versionHistory.rst

@@ -6,6 +6,14 @@
 Version History
 ##################
 
+.. _lsst.ts.wep-8.0.3:
+
+-------------
+8.0.3
+-------------
+
+* Attach locally linear WCSs to DonutStamps.
+
 .. _lsst.ts.wep-8.0.2:
 
 -------------

python/lsst/ts/wep/task/cutOutDonutsBase.py

@@ -29,8 +29,10 @@
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
 import numpy as np
+from lsst.afw.geom import makeSkyWcs
 from lsst.daf.base import PropertyList
 from lsst.fgcmcal.utilities import lookupStaticCalibrations
+from lsst.geom import Point2D, degrees
 from lsst.pipe.base import connectionTypes
 from lsst.ts.wep.cwfs.donutTemplateFactory import DonutTemplateFactory
 from lsst.ts.wep.cwfs.instrument import Instrument
@@ -479,21 +481,33 @@ def cutOutStamps(self, exposure, donutCatalog, defocalType, cameraName):
             else:
                 blendCentroidPositions = np.array([["nan"], ["nan"]], dtype=float).T
 
+            # Get the local linear WCS for the donut stamp
+            # Be careful to get the cd matrix from the linearized WCS instead
+            # of the one from the full WCS.
+            wcs = exposure.wcs
+            centroid_position = Point2D(finalDonutX, finalDonutY)
+            linearTransform = wcs.linearizePixelToSky(centroid_position, degrees)
+            cdMatrix = linearTransform.getLinear().getMatrix()
+            linear_wcs = makeSkyWcs(
+                centroid_position, wcs.pixelToSky(centroid_position), cdMatrix
+            )
+
             donutStamp = DonutStamp(
                 stamp_im=finalStamp,
                 sky_position=lsst.geom.SpherePoint(
                     donutRow["coord_ra"],
                     donutRow["coord_dec"],
                     lsst.geom.radians,
                 ),
-                centroid_position=lsst.geom.Point2D(finalDonutX, finalDonutY),
+                centroid_position=centroid_position,
                 blend_centroid_positions=blendCentroidPositions,
                 detector_name=detectorName,
                 cam_name=cameraName,
                 defocal_type=defocalType.value,
                 # Save defocal offset in mm.
                 defocal_distance=self.instParams["offset"],
                 bandpass=bandpass,
+                archive_element=linear_wcs,
             )
             boundaryT = 1
             maskScalingFactorLocal = 1
@@ -556,4 +570,4 @@ def cutOutStamps(self, exposure, donutCatalog, defocalType, cameraName):
         stampsMetadata["X0"] = np.array(xCornerList)
         stampsMetadata["Y0"] = np.array(yCornerList)
 
-        return DonutStamps(finalStamps, metadata=stampsMetadata)
+        return DonutStamps(finalStamps, metadata=stampsMetadata, use_archive=True)

python/lsst/ts/wep/task/cutOutDonutsCwfsTask.py

@@ -165,8 +165,8 @@ def run(
         extraCatalog, intraCatalog = donutCatalogs
 
         # Get the donut stamps from extra and intra focal images
-        donutStampsExtra = DonutStamps([])
-        donutStampsIntra = DonutStamps([])
+        donutStampsExtra = DonutStamps([], use_archive=True)
+        donutStampsIntra = DonutStamps([], use_archive=True)
 
         for exposure in exposures:
             focusZ = exposure.visitInfo.focusZ

python/lsst/ts/wep/task/donutStamp.py

@@ -281,3 +281,15 @@ def makeMasks(self, inst, model, boundaryT, maskScalingFactorLocal):
         self.mask_pupil = afwImage.Mask(
             np.array(self.comp_im.mask_pupil, dtype=np.int32)
         )
+
+    def getLinearWCS(self):
+        """
+        Get the linear WCS for the stamp.
+
+        Returns
+        -------
+        `lsst.afw.geom.SkyWcs`
+            Linear WCS for the stamp.
+        """
+
+        return self.archive_element

tests/task/test_cutOutDonutsBase.py

@@ -269,6 +269,23 @@ def testCutOutStamps(self):
             donutStamps[0].stamp_im.image.array, expCutOut
         )
 
+        # Check that local linear WCS in archive element is consistent with the
+        # original exposure WCS.
+        exposure_wcs = exposure.wcs
+        for stamp in donutStamps:
+            stamp_wcs = stamp.getLinearWCS()
+            pt = stamp.centroid_position
+            for offset in [(0, 0), (10, -20), (-30, 40), (50, 60)]:
+                pt_offset = pt + lsst.geom.Extent2D(*offset)
+                self.assertFloatsAlmostEqual(
+                    exposure_wcs.pixelToSky(pt_offset)
+                    .separation(stamp_wcs.pixelToSky(pt_offset))
+                    .asArcseconds(),
+                    0.0,
+                    rtol=0.0,
+                    atol=10e-6,  # 10 microarcsecond accurate over stamp region
+                )
+
     def testCutOutStampsBlended(self):
         exposure = self.butler.get(
             "postISRCCD", dataId=self.dataIdExtra, collections=[self.runName]

tests/task/test_donutStamp.py

@@ -180,6 +180,15 @@ def testGetCamera(self):
         with self.assertRaises(ValueError, msg=errMessage):
             donutStamp.getCamera()
 
+    def testGetLinearWCS(self):
+        wcs = lsst.afw.geom.makeSkyWcs(
+            lsst.geom.Point2D(0.0, 0.0),
+            lsst.geom.SpherePoint(0.0, 0.0, lsst.geom.degrees),
+            np.eye(2),
+        )
+        donutStamp = DonutStamp.factory(self.testStamps[0], self.testMetadata, 0, wcs)
+        self.assertEqual(wcs, donutStamp.getLinearWCS())
+
     def testCalcFieldXY(self):
         donutStamp = DonutStamp(
             self.testStamps[0],