Add local linear WCSs to DonutStamps.

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,6 +481,17 @@ 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
+            wcs = exposure.wcs
+            pt = Point2D(finalDonutX, finalDonutY)
+            spt = wcs.pixelToSky(pt)
+            # Be careful to get the cd matrix from the linearized WCS instead
+            # of the one from the full WCS.
+            lin = wcs.linearizePixelToSky(pt, degrees)
+            cd = lin.getLinear().getMatrix()
+            linear_wcs = makeSkyWcs(pt, spt, cd)
+            archive_element = linear_wcs
+
             donutStamp = DonutStamp(
                 stamp_im=finalStamp,
                 sky_position=lsst.geom.SpherePoint(
@@ -494,6 +507,7 @@ def cutOutStamps(self, exposure, donutCatalog, defocalType, cameraName):
                 # Save defocal offset in mm.
                 defocal_distance=self.instParams["offset"],
                 bandpass=bandpass,
+                archive_element=archive_element,
             )
             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

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.archive_element
+            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]