Merge PR#305 and PR#306 (fix generateDonutDirectDetect)

.gitignore

@@ -12,6 +12,7 @@ config.log
 # Built by sconsUtils
 version.py
 /bin/
+python/*.dist-info/
 
 # Pytest
 tests/.tests

doc/versionHistory.rst

@@ -6,6 +6,22 @@
 Version History
 ##################
 
+.. _lsst.ts.wep-13.3.2:
+
+-------------
+13.3.2
+-------------
+
+* Fix generateDonutDirectDetect for null donut selection.
+
+.. _lsst.ts.wep-13.3.1:
+
+-------------
+13.3.1
+-------------
+
+* Add isr configs back into default pipelines.
+
 .. _lsst.ts.wep-13.3.0:
 
 -------------

pipelines/_ingredients/wepDirectDetectScienceGroupPipeline.yaml

@@ -1,6 +1,19 @@
 description: wep direct detect pipeline
 
 tasks:
+  isr:
+    class: lsst.ip.isr.IsrTaskLSST
+    config:
+      # Although we don't have to apply the amp offset corrections, we do want
+      # to compute them for analyzeAmpOffsetMetadata to report on as metrics.
+      doAmpOffset: true
+      ampOffset.doApplyAmpOffset: false
+      # Turn off slow steps in ISR
+      doBrighterFatter: false
+      # Mask saturated pixels,
+      # but turn off quadratic crosstalk because it's currently broken
+      doSaturation: True
+      crosstalk.doQuadraticCrosstalkCorrection: False
   generateDonutDirectDetectTask:
     class: lsst.ts.wep.task.generateDonutDirectDetectTask.GenerateDonutDirectDetectTask
     config:

pipelines/production/comCamSimRapidAnalysisPipeline.yaml

@@ -3,6 +3,9 @@ instrument: lsst.obs.lsst.LsstComCamSim
 imports:
   - $TS_WEP_DIR/pipelines/_ingredients/wepDirectDetectScienceGroupPipeline.yaml
   - $TS_WEP_DIR/pipelines/_ingredients/donutVizGroupPipeline.yaml
+tasks:
+  aggregateDonutTablesTask:
+      class: lsst.donut.viz.AggregateDonutTablesTask
 
 # Define pipeline steps
 subsets:

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

@@ -608,38 +608,39 @@ def cutOutStamps(self, exposure, donutCatalog, defocalType, cameraName):
             # Save MaskedImage to stamp
             finalStamp = finalCutout.getMaskedImage()
 
-            # Save centroid positions as str so we can store in header
-            blendStrX = ""
-            blendStrY = ""
-
-            for blend_cx, blend_cy in zip(
-                catalogMeta["blend_centroid_x"][idx],
-                catalogMeta["blend_centroid_y"][idx],
-            ):
-                blend_final_x = blend_cx + donutRow["xShift"]
-                blend_final_y = blend_cy + donutRow["yShift"]
-                blendStrX += f"{blend_final_x:.2f},"
-                blendStrY += f"{blend_final_y:.2f},"
-            # Remove comma from last entry
-            if len(blendStrX) > 0:
-                blendStrX = blendStrX[:-1]
-                blendStrY = blendStrY[:-1]
-            else:
-                blendStrX = None
-                blendStrY = None
-            finalBlendXList.append(blendStrX)
-            finalBlendYList.append(blendStrY)
-
-            # Prepare blend centroid position information
-            if len(catalogMeta["blend_centroid_x"][idx]) > 0:
-                blendCentroidPositions = np.array(
-                    [
-                        catalogMeta["blend_centroid_x"][idx] + donutRow["xShift"],
-                        catalogMeta["blend_centroid_y"][idx] + donutRow["yShift"],
-                    ]
-                ).T
-            else:
-                blendCentroidPositions = np.array([["nan"], ["nan"]], dtype=float).T
+            # Set that as default, unless overwritten
+            blendCentroidPositions = np.array([["nan"], ["nan"]], dtype=float).T
+
+            if len(catalogMeta["blend_centroid_x"]) > 0:
+                # Save centroid positions as str so we can store in header
+                blendStrX = ""
+                blendStrY = ""
+                for blend_cx, blend_cy in zip(
+                    catalogMeta["blend_centroid_x"][idx],
+                    catalogMeta["blend_centroid_y"][idx],
+                ):
+                    blend_final_x = blend_cx + donutRow["xShift"]
+                    blend_final_y = blend_cy + donutRow["yShift"]
+                    blendStrX += f"{blend_final_x:.2f},"
+                    blendStrY += f"{blend_final_y:.2f},"
+                # Remove comma from last entry
+                if len(blendStrX) > 0:
+                    blendStrX = blendStrX[:-1]
+                    blendStrY = blendStrY[:-1]
+                else:
+                    blendStrX = None
+                    blendStrY = None
+                finalBlendXList.append(blendStrX)
+                finalBlendYList.append(blendStrY)
+
+                # Prepare blend centroid position information
+                if len(catalogMeta["blend_centroid_x"][idx]) > 0:
+                    blendCentroidPositions = np.array(
+                        [
+                            catalogMeta["blend_centroid_x"][idx] + donutRow["xShift"],
+                            catalogMeta["blend_centroid_y"][idx] + donutRow["yShift"],
+                        ]
+                    ).T
 
             # Get the local linear WCS for the donut stamp
             # Be careful to get the cd matrix from the linearized WCS instead

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

@@ -199,17 +199,44 @@ def updateDonutCatalog(self, donutCat, exposure):
         ]
         donutCatUpd["source_flux"] = donutCat["source_flux"] * u.nJy
         fluxSort = np.argsort(donutCatUpd["source_flux"])[::-1]
-        donutCatUpd.meta["blend_centroid_x"] = [
-            donutCat.meta["blend_centroid_x"][idx] for idx in fluxSort
-        ]
-        donutCatUpd.meta["blend_centroid_y"] = [
-            donutCat.meta["blend_centroid_y"][idx] for idx in fluxSort
-        ]
+        # It is possible for catalog to have multiple sources
+        # detected, but with source selection turned off
+        # the QTable metadata of `blend_centroid_x` and
+        # `blend_centroid_y` will be empty.
+        if self.config.doDonutSelection:
+            donutCatUpd.meta["blend_centroid_x"] = [
+                donutCat.meta["blend_centroid_x"][idx] for idx in fluxSort
+            ]
+            donutCatUpd.meta["blend_centroid_y"] = [
+                donutCat.meta["blend_centroid_y"][idx] for idx in fluxSort
+            ]
 
         donutCatUpd.sort("source_flux", reverse=True)
 
         return donutCatUpd
 
+    def emptyTable(self):
+        """Return empty donut table if no donuts got
+        detected or selected.
+
+        Returns
+        -------
+        astropy.table.QTable
+            An empty donut table with correct columns.
+        """
+        donutColumns = [
+            "coord_ra",
+            "coord_dec",
+            "centroid_x",
+            "centroid_y",
+            "detector",
+            "source_flux",
+        ]
+        donutTable = QTable(names=donutColumns)
+        donutTable.meta["blend_centroid_x"] = ""
+        donutTable.meta["blend_centroid_y"] = ""
+        return donutTable
+
     @timeMethod
     def run(self, exposure, camera):
         camName = camera.getName()
@@ -255,6 +282,11 @@ def run(self, exposure, camera):
             # Use the aperture flux with a 70 pixel aperture
             donutTable[f"{bandLabel}_flux"] = donutTable["apFlux70"]
 
+            # Set the required columns to be empty, unless
+            # overwritten by donutSelector below
+            donutTable.meta["blend_centroid_x"] = ""
+            donutTable.meta["blend_centroid_y"] = ""
+
             # Run the donut selector task.
             if self.config.doDonutSelection:
                 self.log.info("Running Donut Selector")
@@ -265,35 +297,30 @@ def run(self, exposure, camera):
                 donutCatSelected.meta["blend_centroid_x"] = donutSelection.blendCentersX
                 donutCatSelected.meta["blend_centroid_y"] = donutSelection.blendCentersY
             else:
-                # if donut selector was not run,
-                # set the required columns to be empty
-                donutTable.meta["blend_centroid_x"] = ""
-                donutTable.meta["blend_centroid_y"] = ""
                 donutCatSelected = donutTable
 
             donutCatSelected.rename_column(f"{bandLabel}_flux", "source_flux")
 
-            # update column names and content
-            donutCatUpd = self.updateDonutCatalog(donutCatSelected, exposure)
-            donutCatUpd["detector"] = np.array(
-                [detectorName] * len(donutCatUpd), dtype=str
-            )
+            # If at least one donut got selected, update the column names
+            # and content
+            if len(donutCatSelected) > 0:
+                donutCatUpd = self.updateDonutCatalog(donutCatSelected, exposure)
+                donutCatUpd["detector"] = np.array(
+                    [detectorName] * len(donutCatUpd), dtype=str
+                )
+            # If no donuts got selected, issue a warning and return an empty
+            # donut table
+            else:
+                self.log.warning(
+                    "No sources selected in the exposure. Returning an empty donut catalog."
+                )
+                donutCatUpd = self.emptyTable()
         else:
 
             self.log.warning(
                 "No sources found in the exposure. Returning an empty donut catalog."
             )
-            donutCatalogColumns = [
-                "coord_ra",
-                "coord_dec",
-                "centroid_x",
-                "centroid_y",
-                "detector",
-                "source_flux",
-            ]
-            donutCatUpd = QTable(names=donutCatalogColumns)
-            donutCatUpd.meta["blend_centroid_x"] = ""
-            donutCatUpd.meta["blend_centroid_y"] = ""
+            donutCatUpd = self.emptyTable()
 
         donutCatUpd = addVisitInfoToCatTable(exposure, donutCatUpd)
 

tests/task/test_generateDonutDirectDetectTask.py

@@ -148,6 +148,36 @@ def testUpdateDonutCatalog(self):
         metaKeys = ["blend_centroid_x", "blend_centroid_y"]
         self.assertCountEqual(donutCatUpd.meta.keys(), metaKeys)
 
+        # test that the catalog created without
+        # source selection will work
+        donutCatNoSel = QTable(data={x: y for x, y in zip(names, data)})
+        donutCatNoSel.meta["blend_centroid_x"] = ""
+        donutCatNoSel.meta["blend_centroid_y"] = ""
+        self.task.config.doDonutSelection = False
+
+        # update the donut catalog
+        donutCatUpd = self.task.updateDonutCatalog(donutCatNoSel, testExposure)
+
+        # check that the new columns are present
+        self.assertCountEqual(newColumns, donutCatUpd.columns)
+
+    def testEmptyTable(self):
+
+        testTable = self.task.emptyTable()
+
+        # Test that there are no rows, but all columns are present
+        self.assertEqual(len(testTable), 0)
+
+        expected_columns = [
+            "coord_ra",
+            "coord_dec",
+            "centroid_x",
+            "centroid_y",
+            "detector",
+            "source_flux",
+        ]
+        self.assertCountEqual(testTable.columns, expected_columns)
+
     def testTaskRun(self):
         """
         Test that the task runs interactively.
@@ -193,9 +223,12 @@ def testTaskRun(self):
             "source_flux",
         ]
         self.assertCountEqual(taskOutNoSrc.donutCatalog.columns, expected_columns)
+
+        # Test that all expected metadata keys are present
+        expected_metakeys = ["blend_centroid_x", "blend_centroid_y", "visit_info"]
         self.assertCountEqual(
             taskOutNoSrc.donutCatalog.meta.keys(),
-            ["blend_centroid_x", "blend_centroid_y", "visit_info"],
+            expected_metakeys,
         )
 
         # Run detection with different sources in each exposure
@@ -236,6 +269,54 @@ def testTaskRun(self):
         self.assertLess(diff_x, tolerance)
         self.assertLess(diff_y, tolerance)
 
+        # Test behavior if no sources get selected
+        # This setting will select no donuts
+        # on that exposure
+        self.task.config.donutSelector.maxFieldDist = 0
+
+        # Run the task
+        taskOut_S10_noSources = self.task.run(
+            exposure_S10,
+            camera,
+        )
+
+        # Test that there are no rows, but all columns are present
+        self.assertCountEqual(
+            taskOut_S10_noSources.donutCatalog.columns, expected_columns
+        )
+
+        # Test that all expected metadata keys are present
+        self.assertCountEqual(
+            taskOut_S10_noSources.donutCatalog.meta.keys(),
+            expected_metakeys,
+        )
+
+        # Test the behavior when source selection is turned off
+        self.task.config.doDonutSelection = False
+        taskOut_S11_noSelection = self.task.run(
+            exposure_S11,
+            camera,
+        )
+        # Check that the expected columns are present
+        self.assertCountEqual(
+            taskOut_S11_noSelection.donutCatalog.columns, expected_columns
+        )
+        # Check that the length of catalogs is as expected
+        outputTableNoSel = taskOut_S11_noSelection.donutCatalog
+        self.assertEqual(len(outputTableNoSel), 3)
+
+        # Test against truth the centroid for all sources
+        result_y = np.sort(outputTableNoSel["centroid_y"].value)
+        truth_y = np.sort(np.array([3196, 2196, 398]))
+        diff_y = np.sum(result_y - truth_y)
+
+        result_x = np.sort(outputTableNoSel["centroid_x"].value)
+        truth_x = np.sort(np.array([3812, 2814, 618]))
+        diff_x = np.sum(result_x - truth_x)
+
+        self.assertLess(diff_x, tolerance)
+        self.assertLess(diff_y, tolerance)
+
     def testTaskRunPipeline(self):
         """
         Test that the task runs in a pipeline.

tests/test_pipelines.py

@@ -0,0 +1,45 @@
+# This file is part of ts_wep.
+#
+# Developed for the LSST Telescope and Site Systems.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import unittest
+from glob import glob
+from pathlib import Path
+
+from lsst.pipe.base import Pipeline
+from lsst.utils import getPackageDir
+
+
+class TestPipeline(unittest.TestCase):
+    """Test the pipeline."""
+
+    def testPipeline(self):
+        packageDir = getPackageDir("ts_wep")
+        # only test production pipelines
+        pipelinePattern = Path(packageDir) / "pipelines" / "production"
+        files = glob(pipelinePattern.as_posix() + "/*.yaml")
+        for filename in files:
+            pipeline = Pipeline.fromFile(filename)
+            self.assertIsInstance(pipeline, Pipeline)
+
+
+if __name__ == "__main__":
+    # Do the unit test
+    unittest.main()