Tickets/DM-52459: add donut_id

doc/news/DM-52459.misc.md

@@ -0,0 +1 @@
+Add donut_id generated during the creation of donut table.

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

@@ -299,6 +299,8 @@ def initZkTable(self) -> QTable:
             ("extra_frac_bad_pix", "<f4"),
             ("intra_max_power_grad", "<f4"),
             ("extra_max_power_grad", "<f4"),
+            ("intra_donut_id", "<U21"),
+            ("extra_donut_id", "<U21"),
         ]
         for j in self.nollIndices:
             dtype.append((f"Z{j}", "<f4"))
@@ -362,6 +364,8 @@ def createZkTable(self, zkCoeffRaw: pipeBase.Struct) -> QTable:
                 "extra_frac_bad_pix": np.nan,
                 "intra_max_power_grad": np.nan,
                 "extra_max_power_grad": np.nan,
+                "intra_donut_id": "",
+                "extra_donut_id": "",
             }
         )
         for i, (intra, extra, zk) in enumerate(
@@ -400,15 +404,15 @@ def createZkTable(self, zkCoeffRaw: pipeBase.Struct) -> QTable:
             row["extra_field"] = extraAngle
             row["intra_centroid"] = intraCentroid
             row["extra_centroid"] = extraCentroid
-            for key in ["MAG", "SN", "ENTROPY", "FRAC_BAD_PIX", "MAX_POWER_GRAD"]:
+            for key in ["MAG", "SN", "ENTROPY", "FRAC_BAD_PIX", "MAX_POWER_GRAD", "DONUT_ID"]:
                 for stamps, foc in [
                     (self.stampsIntra, "intra"),
                     (self.stampsExtra, "extra"),
                 ]:
                     if len(stamps) > 0 and key in stamps.metadata:
                         val = stamps.metadata.getArray(key)[i]
                     else:
-                        val = np.nan
+                        val = "" if key == "DONUT_ID" else np.nan
                     row[f"{foc}_{key.lower()}"] = val
             zkTable.add_row(row)
 
@@ -493,6 +497,7 @@ def empty(
             "RADIUS",
             "RADIUS_FAIL_FLAG",
             "DEFOCAL_TYPE",
+            "DONUT_ID",
         ]
         if qualityTable is None:
             donutQualityTable = QTable({name: [] for name in qualityTableCols})

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

@@ -753,6 +753,7 @@ def cutOutStamps(
                     donutRow["coord_dec"].value,
                     lsst.geom.radians,
                 ),
+                donut_id=donutRow["donut_id"],
                 centroid_position=centroid_position,
                 blend_centroid_positions=blendCentroidPositions,
                 detector_name=detectorName,
@@ -800,6 +801,9 @@ def cutOutStamps(
             stampsMetadata["MAG"] = (donutCatalog[fluxLabel].value * u.nJy).to_value(u.ABmag)
         else:
             stampsMetadata["MAG"] = np.array([])
+
+        # Save donut id
+        stampsMetadata["DONUT_ID"] = np.array(donutCatalog["donut_id"].value)
         # Save the original centroid values
         stampsMetadata["CENT_X0"] = np.array(donutCatalog["centroid_x"].value)
         stampsMetadata["CENT_Y0"] = np.array(donutCatalog["centroid_y"].value)

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

@@ -74,6 +74,8 @@ class DonutStamp(AbstractStamp):
         are available camera names currently.
     bandpass : `str`
         The bandpass for the stamp image.
+    donut_id: `str`
+        The donut id defined at detection stage.
     archive_element : `afwTable.io.Persistable`, optional
         Archive element (e.g. Transform or WCS) associated with this stamp.
         (the default is None.)
@@ -83,6 +85,7 @@ class DonutStamp(AbstractStamp):
         camera coordinate system (CCS), with the CWFSs rotated to the same
         orientation as the science sensors. It is this object that will be used
         to interface with the wavefront estimator.
+
     """
 
     stamp_im: afwImage.MaskedImageF
@@ -94,6 +97,7 @@ class DonutStamp(AbstractStamp):
     detector_name: str
     cam_name: str
     bandpass: str
+    donut_id: str
     archive_element: Optional[afwTable.io.Persistable] = None
     wep_im: Image = field(init=False)
 
@@ -178,6 +182,11 @@ def factory(
             # If this is an old version of the stamps without bandpass
             # information then an empty string ("") will be set as default.
             bandpass=(metadata.getArray("BANDPASS")[index] if metadata.get("BANDPASS") is not None else ""),
+            # "DONUT_ID" is the donut identification number
+            # created at detection stage. If this is an old
+            # version of stamps without donut id,
+            # an empty string will be set as default
+            donut_id=(metadata.getArray("DONUT_ID")[index] if metadata.get("DONUT_ID") is not None else ""),
         )
 
     def getCamera(self) -> Camera:

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

@@ -211,6 +211,12 @@ def selectStamps(self, donutStamps: DonutStamps) -> pipeBase.Struct:
         # initiate these by selecting all donuts
         entropySelect = np.ones(len(donutStamps), dtype="bool")
 
+        # Collect donut id values
+        if "DONUT_ID" in list(donutStamps.metadata):
+            donutId = donutStamps.metadata.getArray("DONUT_ID")
+        else:
+            donutId = np.zeros(len(donutStamps), dtype="str")
+
         # Collect the entropy information if available
         entropyValue = np.full(len(donutStamps), np.nan)
         if "ENTROPY" in list(donutStamps.metadata):
@@ -320,6 +326,7 @@ def selectStamps(self, donutStamps: DonutStamps) -> pipeBase.Struct:
                 maxPowerGradSelect,
                 donutRadii,
                 selected,
+                donutId,
             ],
             names=[
                 "SN",
@@ -332,6 +339,7 @@ def selectStamps(self, donutStamps: DonutStamps) -> pipeBase.Struct:
                 "MAX_POWER_GRAD_SELECT",
                 "RADIUS",
                 "FINAL_SELECT",
+                "DONUT_ID",
             ],
         )
         # Add all configuration used for selection criteria as metadata;

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

@@ -59,6 +59,8 @@ def _refresh_metadata(self) -> None:
         self.metadata["DFC_DIST"] = [dfc_dist for dfc_dist in defocal_distances]
         bandpasses = self.getBandpasses()
         self.metadata["BANDPASS"] = [bandpass for bandpass in bandpasses]
+        donut_ids = self.getIds()
+        self.metadata["DONUT_ID"] = [idx for idx in donut_ids]
 
     def getSkyPositions(self) -> list:
         """
@@ -181,6 +183,17 @@ def getBandpasses(self) -> list:
         """
         return [stamp.bandpass for stamp in self]
 
+    def getIds(self) -> list:
+        """
+        Get the id for each stamp.
+
+        Returns
+        -------
+        list [str]
+            Id name for each stamp.
+        """
+        return [stamp.donut_id for stamp in self]
+
     def append(self, newStamp: DonutStamp) -> None:
         """Add an additional stamp.
 

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

@@ -269,4 +269,12 @@ def addVisitInfoToCatTable(exposure: Exposure, donutCat: QTable) -> QTable:
 
     donutCat.meta["visit_info"] = catVisitInfo
 
+    # add donutId once sources are sorted by brightness
+    detId = exposure.detector.getId()
+    visit = exposure.visitInfo.id
+    donutId = [
+        str(visit) + "_" + str(detId).zfill(3) + "_" + str(idx).zfill(3) for idx in np.arange(len(donutCat))
+    ]
+    donutCat["donut_id"] = donutId
+
     return donutCat

tests/task/test_calcZernikesTieTaskCwfs.py

@@ -264,6 +264,10 @@ def testUnevenPairs(self) -> None:
         # Increase length of extra list
         stampsExtra.extend([stampsExtra[0]])
 
+        # Refresh metadata ensures arrays
+        # are same length as new donut stamp list
+        stampsExtra._refresh_metadata()
+
         # Now estimate Zernikes
         self.task.run(stampsExtra, stampsIntra, self.intrinsicTables)
 

tests/task/test_calcZernikesTieTaskScienceSensor.py

@@ -193,6 +193,8 @@ def testCalcZernikes(self) -> None:
         # verify remaining desired columns exist in zernikes table
         desired_colnames = [
             "used",
+            "intra_donut_id",
+            "extra_donut_id",
             "intra_field",
             "extra_field",
             "intra_centroid",
@@ -244,6 +246,7 @@ def testCalcZernikes(self) -> None:
             "MAX_POWER_GRAD_SELECT",
             "FINAL_SELECT",
             "DEFOCAL_TYPE",
+            "DONUT_ID",
             "RADIUS",
         ]
         np.testing.assert_array_equal(np.sort(colnames), np.sort(desired_colnames))

tests/task/test_calcZernikesUnpairedTask.py

@@ -222,6 +222,8 @@ def testTable(self) -> None:
                     "extra_frac_bad_pix",
                     "intra_max_power_grad",
                     "extra_max_power_grad",
+                    "intra_donut_id",
+                    "extra_donut_id",
                 ]
                 self.assertLessEqual(set(desired_colnames), set(structNormal.zernikes.colnames))
 
@@ -260,6 +262,7 @@ def testTable(self) -> None:
                     "FINAL_SELECT",
                     "DEFOCAL_TYPE",
                     "RADIUS",
+                    "DONUT_ID",
                 ]
                 np.testing.assert_array_equal(np.sort(colnames), np.sort(desired_colnames))
 

tests/task/test_cutOutDonutsBase.py

@@ -25,6 +25,7 @@
 from astropy.table import QTable
 from scipy.signal import correlate
 
+import lsst.geom
 import lsst.utils.tests
 from lsst.afw import image as afwImage
 from lsst.daf.base import PropertySet
@@ -434,6 +435,7 @@ def testCutOutStampsTaskRunNormal(self) -> None:
             "BORESIGHT_RA_RAD",
             "BORESIGHT_DEC_RAD",
             "BANDPASS",
+            "DONUT_ID",
         ]
         self.assertCountEqual(metadata, expectedMetadata)
 

tests/task/test_donutStamp.py

@@ -68,6 +68,7 @@ def _makeStamps(self, nStamps: int, stampSize: int, testDefaults: bool = False)
         dfcTypes[:halfStampIdx] = [DefocalType.Intra.value] * halfStampIdx
         dfcDists = np.ones(nStamps) * 1.25
         bandpass = ["r"] * nStamps
+        donutId = np.arange(nStamps)
 
         metadata = PropertyList()
         metadata["RA_DEG"] = ras
@@ -77,6 +78,7 @@ def _makeStamps(self, nStamps: int, stampSize: int, testDefaults: bool = False)
         metadata["DET_NAME"] = detectorNames
         metadata["CAM_NAME"] = camNames
         metadata["DFC_TYPE"] = dfcTypes
+        metadata["DONUT_ID"] = donutId
         if testDefaults is False:
             metadata["DFC_DIST"] = dfcDists
             metadata["BLEND_CX"] = blendCentX
@@ -122,6 +124,9 @@ def testFactory(self) -> None:
             bandpass = donutStamp.bandpass
             self.assertEqual(bandpass, "r")
 
+            donutId = donutStamp.donut_id
+            self.assertEqual(i, donutId)
+
             self.assertIsInstance(donutStamp.wep_im, Image)
 
     def testFactoryMetadataDefaults(self) -> None:
@@ -145,6 +150,9 @@ def testFactoryMetadataDefaults(self) -> None:
             # Test default bandpass value is empty string
             bandpass = donutStamp.bandpass
             self.assertEqual(bandpass, "")
+            # Test default donutId value is empty string
+            donutId = donutStamp.donut_id
+            self.assertEqual(donutId, i)
 
     def testGetCamera(self) -> None:
         donutStamp = DonutStamp.factory(self.testStamps[0], self.testMetadata, 0)
@@ -189,6 +197,7 @@ def testCalcFieldXY(self) -> None:
             "R22_S11",
             "LSSTCam",
             "r",
+            "",
         )
         np.testing.assert_array_almost_equal(donutStamp.calcFieldXY(), (0, 0))
 
@@ -213,6 +222,7 @@ def testCalcFieldXY(self) -> None:
                     detName,
                     "LSSTCam",
                     "r",
+                    "",
                 )
                 fieldAngle = donutStamp.calcFieldXY()
                 self.assertEqual(fieldAngle[0], np.degrees(trueFieldAngleX))
@@ -229,6 +239,7 @@ def testMakeMask(self) -> None:
             "R22_S11",
             "LSSTCam",
             "r",
+            "",
         )
 
         # Check that mask is empty at start
@@ -284,6 +295,7 @@ def _testWepImage(
                 detName,
                 "LSSTCam",
                 "r",
+                "",
             )
 
             # Make the mask

tests/task/test_donutStampSelectorTask.py

@@ -190,12 +190,12 @@ def testSelectStamps(self) -> None:
         # test custom SNR thresholds
         self.config.selectWithEntropy = False
         self.config.useCustomSnLimit = True
-        minSignalToNoise = 1000.0
+        minSignalToNoise = 1700.0
         self.config.minSignalToNoise = minSignalToNoise
         task = DonutStampSelectorTask(config=self.config, name="SN Task")
         selection = task.selectStamps(donutStampsIntra)
         donutsQuality = selection.donutsQuality
-        self.assertEqual(np.sum(donutsQuality["SN_SELECT"]), 3)
+        self.assertEqual(np.sum(donutsQuality["SN_SELECT"]), 0)
 
         # test that the SN of selected donuts is indeed above the threshold
         for v in donutsQuality["SN"][donutsQuality["SN_SELECT"]]:
@@ -230,6 +230,11 @@ def testSelectStamps(self) -> None:
         selection = task.selectStamps(donutStampsIntra)
         self.assertEqual(np.sum(selection.donutsQuality["FINAL_SELECT"]), 3)
 
+        # test whether DonutStamps metadata get copied over if present
+        for metakey in ["DONUT_ID", "RADIUS"]:
+            if metakey in list(donutStampsIntra.metadata):
+                self.assertIn(metakey, selection.donutsQuality.columns)
+
     def testTaskRun(self) -> None:
         donutStampsIntra = self.butler.get(
             "donutStampsIntra", dataId=self.dataIdExtra, collections=[self.baseRunName]

tests/task/test_donutStamps.py

@@ -68,6 +68,7 @@ def _makeDonutStamps(
         dfcTypes[:halfStampIdx] = [DefocalType.Intra.value] * halfStampIdx
         dfcDists = np.ones(nStamps) * 1.5
         bandpass = ["r"] * nStamps
+        donutId = ["0000"] * nStamps
 
         # Test mixture of donuts with blends and those without
         blendCentX[-1] = "nan"
@@ -85,6 +86,7 @@ def _makeDonutStamps(
         metadata["DFC_TYPE"] = dfcTypes
         metadata["DFC_DIST"] = dfcDists
         metadata["BANDPASS"] = bandpass
+        metadata["DONUT_ID"] = donutId
 
         donutStampList = [DonutStamp.factory(stampList[idx], metadata, idx) for idx in range(nStamps)]
 
@@ -115,6 +117,7 @@ def _roundtrip(self, donutStamps: DonutStamps) -> None:
                 self.assertEqual(stamp1.defocal_type, stamp2.defocal_type)
                 self.assertEqual(stamp1.defocal_distance, stamp2.defocal_distance)
                 self.assertEqual(stamp1.bandpass, stamp2.bandpass)
+                self.assertEqual(stamp1.donut_id, stamp2.donut_id)
 
     def testGetSkyPositions(self) -> None:
         skyPos = self.donutStamps.getSkyPositions()
@@ -171,6 +174,10 @@ def testGetBandpass(self) -> None:
         bandpasses = self.donutStamps.getBandpasses()
         self.assertListEqual(bandpasses, ["r"] * self.nStamps)
 
+    def testGetDonutId(self) -> None:
+        donutId = self.donutStamps.getIds()
+        self.assertListEqual(donutId, ["0000"] * self.nStamps)
+
     def testAppend(self) -> None:
         """Test ability to append to a Stamps object"""
         self.donutStamps.append(self.donutStamps[0])
@@ -180,7 +187,7 @@ def testAppend(self) -> None:
         self._roundtrip(self.donutStamps)
         # check if appending something other than a DonutStamp raises
         with self.assertRaises(ValueError) as context:
-            self.donutStamps.append("hello world")
+            self.donutStamps.append("hello world")  # type: ignore[arg-type]
         self.assertEqual("Objects added must be a DonutStamp object.", str(context.exception))
 
     def testExtend(self) -> None:

tests/task/test_generateDonutCatalogUtils.py

@@ -367,3 +367,20 @@ def testAddVisitInfoToCatTable(self) -> None:
         self.assertTrue(isinstance(catTableWithMeta.meta["visit_info"], dict))
         # Test columns are all present
         self.assertCountEqual(visitInfoKeys, catTableWithMeta.meta["visit_info"].keys())
+
+        # Test that donut_id column is present
+        self.assertIn("donut_id", catTableWithMeta.colnames)
+
+        # Test that donut_id has correct length
+        self.assertEqual(len(catTableWithMeta["donut_id"]), len(fieldObjects))
+
+        # Test that donut_id has correct format: {visit}_{detId:03d}_{idx:03d}
+        visit = testExposure.visitInfo.id
+        detId = testExposure.detector.getId()
+
+        for idx, donut_id in enumerate(catTableWithMeta["donut_id"]):
+            expected_id = f"{visit}_{str(detId).zfill(3)}_{str(idx).zfill(3)}"
+            self.assertEqual(donut_id, expected_id)
+
+        # Test that all donut_ids are unique
+        self.assertEqual(len(set(catTableWithMeta["donut_id"])), len(catTableWithMeta))