Adding a basic source catalog from Source Detection Step

CHANGES.rst

@@ -91,7 +91,9 @@ source_detection
 
 - Support for PSF fitting (optional) for accurate centroids. [#841, #984]
 
-- Save source catalog to a structured array. [#987]
+- Save tweakreg source catalog to an astropy Table. [#987, #1029]
+
+- Save basic source catalog with morphology and PSF fit metrics. [#1029]
 
 stpipe
 ------

romancal/lib/basic_utils.py

@@ -103,24 +103,3 @@ def __exit__(self, et, ev, tb):
         if self.handler and self.close:
             self.handler.close()
         # implicit return of None => don't swallow exceptions
-
-
-def recarray_to_ndarray(x, to_dtype="<f8"):
-    """
-    Convert a structured array to a 2D ndarray.
-
-    Parameters
-    ----------
-    x : np.record
-        Structured array
-    to_dtype : str
-        Cast all columns in `x` to this dtype in the output ndarray.
-
-    Returns
-    -------
-    array : np.ndarray
-        Numpy array (without labeled columns).
-    """
-    names = x.dtype.names
-    astype = [(name, to_dtype) for name in names]
-    return np.asarray(x.astype(astype).view(to_dtype).reshape((-1, len(names))))

romancal/lib/tests/test_basic_utils.py

@@ -1,10 +1,8 @@
 """Test basic utils"""
 
-import numpy as np
 import pytest
-from astropy.table import Table
 
-from romancal.lib.basic_utils import bytes2human, recarray_to_ndarray
+from romancal.lib.basic_utils import bytes2human
 
 test_data = [
     (1000, "1000B"),
@@ -19,21 +17,3 @@ def test_bytes2human(input_data, result):
     """test the basic conversion"""
 
     assert bytes2human(input_data) == result
-
-
-def test_structured_array_utils():
-    arrays = [np.arange(0, 10), np.arange(10, 20), np.arange(30, 40)]
-    names = "a, b, c"
-
-    recarr0 = np.core.records.fromarrays(
-        arrays, names=names, formats=[arr.dtype for arr in arrays]
-    )
-    round_tripped = recarray_to_ndarray(recarr0)
-    assert np.all(round_tripped == np.column_stack(arrays).astype("<f8"))
-
-    astropy_table = Table(recarr0)
-    assert np.all(
-        np.array(arrays) == np.array([col.data for col in astropy_table.itercols()])
-    )
-
-    assert np.all(astropy_table.as_array() == recarr0)

romancal/source_detection/source_detection_step.py

@@ -9,7 +9,7 @@
 import numpy as np
 from asdf import AsdfFile
 from astropy.stats import SigmaClip
-from astropy.table import Table
+from astropy.table import Table, join
 from photutils.background import (
     Background2D,
     MeanBackground,
@@ -150,12 +150,12 @@
             columns = ["id", "xcentroid", "ycentroid", "flux"]
 
             if sources:
-                catalog = sources[columns]
-                log.info(f"Found {len(catalog)} sources.")
+                dao_catalog = sources[columns]
+                log.info(f"Found {len(dao_catalog)} sources.")
             else:
                 # if no sources were detected, return an empty table
                 self.log.warning("No sources detected, returning empty catalog.")
-                catalog = Table(
+                dao_catalog = Table(
                     names=columns,
                     dtype=(int, np.float64, np.float64, np.float64),
                 )
@@ -178,13 +178,11 @@
                 psf_photometry_table, photometry = psf.fit_psf_to_image_model(
                     image_model=input_model,
                     psf_model=gridded_psf_model,
-                    x_init=catalog["xcentroid"],
-                    y_init=catalog["ycentroid"],
+                    x_init=dao_catalog["xcentroid"],
+                    y_init=dao_catalog["ycentroid"],
                     exclude_out_of_bounds=True,
                 )
 
-            catalog_as_recarray = catalog.as_array()
-
             # create meta.source detection section in file
             # if save_catalogs is True, this will be updated with the
             # attribute 'tweakreg_catalog_name' to point to the location
@@ -204,11 +202,11 @@
                 log.info(f"Saving catalog to file: {cat_filename}.")
 
                 if self.output_cat_filetype == "asdf":
-                    tree = {"tweakreg_catalog": catalog_as_recarray}
+                    tree = {"tweakreg_catalog": dao_catalog}
                     ff = AsdfFile(tree)
                     ff.write_to(cat_filename)
                 else:
-                    catalog.write(cat_filename, format="ascii.ecsv", overwrite=True)
+                    dao_catalog.write(cat_filename, format="ascii.ecsv", overwrite=True)
 
                 input_model.meta.source_detection[
                     "tweakreg_catalog_name"
@@ -218,22 +216,33 @@
                     # PSF photometry centroid results are stored in an astropy table
                     # in columns "x_fit" and "y_fit", which we'll rename for
                     # compatibility with tweakreg:
-                    catalog = psf_photometry_table.copy()
-
-                    # rename the PSF photometry table's columns to match the
-                    # expectated columns in tweakreg:
+                    source_catalog = join(
+                        psf_photometry_table,
+                        sources,
+                        keys_left=[
+                            "x_init",
+                            "y_init",
+                        ],  # match these inputs for PSF fitting...
+                        keys_right=[
+                            "xcentroid",
+                            "ycentroid",
+                        ],  # ...with centroid outputs from DAOStarFinder
+                        table_names=["psf", "dao"],
+                    )
+
+                    # DAOStarFinder returns columns "xcentroid" and "ycentroid", which
+                    # we'll append with "_dao" to make their meaning clear
                     for old_name, new_name in zip(
-                        ["x_fit", "y_fit"], ["xcentroid", "ycentroid"]
+                        ["xcentroid", "ycentroid", "flux"],
+                        ["xcentroid_dao", "ycentroid_dao", "flux_dao"],
                     ):
-                        catalog.rename_column(old_name, new_name)
+                        source_catalog.rename_column(old_name, new_name)
 
-                    input_model.meta.source_detection["psf_catalog"] = catalog
+                    input_model.meta.source_detection["source_catalog"] = source_catalog
                 else:
                     # only attach catalog to file if its being passed to the next step
                     # and save_catalogs is false, since it is not in the schema
-                    input_model.meta.source_detection[
-                        "tweakreg_catalog"
-                    ] = catalog_as_recarray
+                    input_model.meta.source_detection["tweakreg_catalog"] = dao_catalog
             input_model.meta.cal_step["source_detection"] = "COMPLETE"
 
         # just pass input model to next step - catalog is stored in meta

romancal/source_detection/tests/test_source_detection_step.py

@@ -1,14 +1,13 @@
 """
- Unit tests for the Roman source detection step code
+Unit tests for the Roman source detection step code.
 """
 
 from copy import deepcopy
 
+import astropy.units as u
 import numpy as np
 import pytest
-from astropy import units as u
 
-from romancal.lib.basic_utils import recarray_to_ndarray
 from romancal.lib.tests.test_psf import add_sources, setup_inputs
 from romancal.source_detection import SourceDetectionStep
 
@@ -39,19 +38,19 @@ def test_dao_vs_psf(self, seed=0):
         source_detect.scalar_threshold = 100
         source_detect.peakmax = None
         dao_result = source_detect.process(image_model)
-        idx, x_dao, y_dao, amp_dao = recarray_to_ndarray(
-            dao_result.meta.source_detection.tweakreg_catalog
-        ).T
+        idx, x_dao, y_dao, amp_dao = list(
+            dao_result.meta.source_detection.tweakreg_catalog.itercols()
+        )
 
         # check that all injected targets are found by DAO:
         assert len(x_dao) == len(x_true)
 
         source_detect.fit_psf = True
         psf_result = source_detect.process(image_model)
-        psf_catalog = psf_result.meta.source_detection.psf_catalog
+        source_catalog = psf_result.meta.source_detection.source_catalog
 
-        extract_columns = ["xcentroid", "x_err", "ycentroid", "y_err", "flux_fit"]
-        x_psf, x_err, y_psf, y_err, amp_psf = psf_catalog[extract_columns].itercols()
+        extract_columns = ["x_fit", "x_err", "y_fit", "y_err", "flux_fit"]
+        x_psf, x_err, y_psf, y_err, amp_psf = source_catalog[extract_columns].itercols()
 
         # check that typical PSF centroids are more accurate than DAO centroids:
         assert np.median(np.abs(x_dao - x_true)) > np.median(np.abs(x_psf - x_true))
@@ -64,3 +63,16 @@ def test_dao_vs_psf(self, seed=0):
 
         # check that typical residuals are consistent with their errors:
         assert np.median(np.abs(x_psf - x_true) / x_err) < 2
+
+        # check that the source catalog contains some expected columns:
+        expected_cols = {
+            # PSF fitting outputs:
+            "x_init",
+            "x_fit",
+            "flux_fit",
+            "qfit",
+            # DAOStarFinder outputs:
+            "xcentroid_dao",
+            "sharpness",
+        }
+        assert expected_cols.issubset(source_catalog.colnames)

romancal/tweakreg/tweakreg_step.py

@@ -174,10 +174,9 @@ def process(self, input):
                     image_model.meta.source_detection, "tweakreg_catalog_name"
                 )
                 if is_tweakreg_catalog_present:
-                    # read catalog from structured array
-                    catalog = Table(
-                        np.asarray(image_model.meta.source_detection.tweakreg_catalog)
-                    )
+                    # read catalog produced in the Source Detection Step, ensure that
+                    # it is an astropy Table:
+                    catalog = Table(image_model.meta.source_detection.tweakreg_catalog)
                 elif is_tweakreg_catalog_name_present:
                     catalog = Table.read(
                         image_model.meta.source_detection.tweakreg_catalog_name,
@@ -244,7 +243,7 @@ def process(self, input):
                     )
 
             # set meta.tweakreg_catalog
-            image_model.meta["tweakreg_catalog"] = catalog.as_array()
+            image_model.meta["tweakreg_catalog"] = catalog
 
             nsources = len(catalog)
             if nsources == 0:
@@ -556,7 +555,7 @@ def _imodel2wcsim(self, image_model):
                 # a string with the name of the catalog was provided
                 catalog = Table.read(catalog, format=catalog_format)
             else:
-                # catalog is a structured array, convert to astropy table:
+                # catalog is an astropy Table or structured array, convert to astropy table:
                 catalog = Table(catalog)
 
             catalog.meta["name"] = (