Copy input datamodel in source_catalog step (#1457)

changes/1457.source_catalog.rst

@@ -0,0 +1,2 @@
+The data and err array of the input datamodel to the source_catalog step
+are now copied so that they are left completely unchanged by the step.

romancal/source_catalog/source_catalog.py

@@ -175,19 +175,6 @@ def convert_l2_to_sb(self):
         self.model["err"] *= self.l2_conv_factor
         self.convolved_data *= self.l2_conv_factor
 
-    def convert_sb_to_l2(self):
-        """
-        Convert the data and error arrays from MJy/sr (surface
-        brightness) to DN/s (level-2 units).
-
-        This is the inverse operation of `convert_l2_to_sb`.
-        """
-        # the conversion in done in-place to avoid making copies of the data;
-        # use a dictionary to set the value to avoid on-the-fly validation
-        self.model["data"] /= self.l2_conv_factor
-        self.model["err"] /= self.l2_conv_factor
-        self.convolved_data /= self.l2_conv_factor
-
     def convert_sb_to_flux_density(self):
         """
         Convert the data and error arrays from MJy/sr (surface
@@ -206,22 +193,6 @@ def convert_sb_to_flux_density(self):
         self.model["err"] <<= self.sb_to_flux.unit
         self.convolved_data <<= self.sb_to_flux.unit
 
-    def convert_flux_density_to_sb(self):
-        """
-        Convert the data and error arrays from flux density units to
-        MJy/sr (surface brightness).
-
-        The units are also removed from the arrays.
-
-        This is the inverse operation of `convert_sb_to_flux_density`.
-        """
-        self.model["data"] /= self.sb_to_flux
-        self.model["err"] /= self.sb_to_flux
-        self.convolved_data /= self.sb_to_flux
-        self.model["data"] = self.model["data"].value
-        self.model["err"] = self.model["err"].value
-        self.convolved_data = self.convolved_data.value
-
     def convert_flux_to_abmag(self, flux, flux_err):
         """
         Convert flux (and error) to AB magnitude (and error).
@@ -1199,11 +1170,4 @@ def catalog(self):
         # split SkyCoord columns into separate RA and Dec columns
         catalog = self._split_skycoord(catalog)
 
-        # restore units on input model back to MJy/sr
-        self.convert_flux_density_to_sb()
-
-        # restore L2 data units back to DN/s
-        if isinstance(self.model, ImageModel):
-            self.convert_sb_to_l2()
-
         return catalog

romancal/source_catalog/source_catalog_step.py

@@ -6,6 +6,7 @@
 from astropy.table import Table
 from roman_datamodels import datamodels, maker_utils
 from roman_datamodels.datamodels import ImageModel, MosaicModel
+from roman_datamodels.maker_utils import mk_datamodel
 
 from romancal.source_catalog.background import RomanBackground
 from romancal.source_catalog.detection import convolve_data, make_segmentation_image
@@ -45,15 +46,39 @@ class SourceCatalogStep(RomanStep):
         fit_psf = boolean(default=True)      # fit source PSFs for accurate astrometry?
     """
 
-    def process(self, input_model):
-        if isinstance(input_model, datamodels.DataModel):
-            model = input_model
+    def process(self, step_input):
+        if isinstance(step_input, datamodels.DataModel):
+            input_model = step_input
         else:
-            model = datamodels.open(input_model)
+            input_model = datamodels.open(step_input)
 
-        if not isinstance(model, (ImageModel, MosaicModel)):
+        if not isinstance(input_model, (ImageModel, MosaicModel)):
             raise ValueError("The input model must be an ImageModel or MosaicModel.")
 
+        # Copy the data and error arrays to avoid modifying the input
+        # model. We use mk_datamodel to copy *only* the data and err
+        # arrays. The metadata and dq and weight arrays are not copied
+        # because they are not modified in this step. The other model
+        # arrays (e.g., var_rnoise) are not currently used by this step.
+        if isinstance(input_model, ImageModel):
+            model = mk_datamodel(
+                ImageModel,
+                meta=input_model.meta,
+                shape=(0, 0),
+                data=input_model.data.copy(),
+                err=input_model.err.copy(),
+                dq=input_model.dq,
+            )
+        elif isinstance(input_model, MosaicModel):
+            model = mk_datamodel(
+                MosaicModel,
+                meta=input_model.meta,
+                shape=(0, 0),
+                data=input_model.data.copy(),
+                err=input_model.err.copy(),
+                weight=input_model.weight,
+            )
+
         if isinstance(model, ImageModel):
             cat_model = datamodels.SourceCatalogModel
         else:
@@ -111,40 +136,34 @@ def process(self, input_model):
             # put the resulting catalog in the model
             source_catalog_model.source_catalog = catobj.catalog
 
-        # add back background to data so input model is unchanged
-        # (in case of interactive use)
-        model.data += bkg.background
-
-        # create catalog filename
-        # N.B.: self.save_model will determine whether to use fully qualified path or not
-        output_catalog_name = self.make_output_path(
-            basepath=model.meta.filename, suffix="cat"
-        )
-
-        if isinstance(model, ImageModel):
-            update_metadata(model, output_catalog_name)
-
-        output_model = model
-
-        # always save segmentation image + catalog
+        # always save the segmentation image and source catalog
         self.save_base_results(segment_img, source_catalog_model)
 
-        # return the updated model or the source catalog object
+        # Return the source catalog object or the input model. If the
+        # input model is an ImageModel, the metadata is updated with the
+        # source catalog filename.
         if getattr(self, "return_updated_model", False):
-            # setting the suffix to something else to prevent
-            # step from  overwriting source catalog file with
-            # a datamodel
+            # define the catalog filename; self.save_model will
+            # determine whether to use a fully qualified path
+            output_catalog_name = self.make_output_path(
+                basepath=model.meta.filename, suffix="cat"
+            )
+
+            # set the suffix to something else to prevent the step from
+            # overwriting the source catalog file with a datamodel
             self.suffix = "sourcecatalog"
-            # return DataModel
-            result = output_model
+
+            if isinstance(input_model, ImageModel):
+                update_metadata(input_model, output_catalog_name)
+
+            result = input_model
         else:
-            # return SourceCatalogModel
             result = source_catalog_model
 
         return result
 
     def save_base_results(self, segment_img, source_catalog_model):
-        # save the segmentation map and
+        # save the segmentation map and source catalog
         output_filename = (
             self.output_file
             if self.output_file is not None
@@ -168,6 +187,7 @@ def save_base_results(self, segment_img, source_catalog_model):
                 suffix="segm",
                 force=True,
             )
+
         # save the source catalog
         self.save_model(
             source_catalog_model, output_file=output_filename, suffix="cat", force=True

romancal/source_catalog/tests/test_source_catalog.py

@@ -6,7 +6,7 @@
 import pytest
 from astropy.modeling.models import Gaussian2D
 from astropy.table import Table
-from numpy.testing import assert_allclose
+from numpy.testing import assert_equal
 from photutils.segmentation import SegmentationImage
 from roman_datamodels import datamodels as rdm
 from roman_datamodels.datamodels import (
@@ -300,8 +300,8 @@ def test_l2_input_model_unchanged(image_model, tmp_path):
         save_results=False,
     )
 
-    assert_allclose(original_data, image_model.data, atol=5.0e-5)
-    assert_allclose(original_err, image_model.err, atol=5.0e-5)
+    assert_equal(original_data, image_model.data)
+    assert_equal(original_err, image_model.err)
 
 
 @pytest.mark.webbpsf
@@ -324,8 +324,8 @@ def test_l3_input_model_unchanged(mosaic_model, tmp_path):
         save_results=False,
     )
 
-    assert_allclose(original_data, mosaic_model.data, atol=5.0e-5)
-    assert_allclose(original_err, mosaic_model.err, atol=5.0e-5)
+    assert_equal(original_data, mosaic_model.data)
+    assert_equal(original_err, mosaic_model.err)
 
 
 @pytest.mark.webbpsf