diff --git a/changes/1457.source_catalog.rst b/changes/1457.source_catalog.rst new file mode 100644 index 000000000..16581605f --- /dev/null +++ b/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. diff --git a/romancal/source_catalog/source_catalog.py b/romancal/source_catalog/source_catalog.py index 6357d4831..b22a5fe07 100644 --- a/romancal/source_catalog/source_catalog.py +++ b/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 diff --git a/romancal/source_catalog/source_catalog_step.py b/romancal/source_catalog/source_catalog_step.py index 9946d037b..e1a230b2f 100644 --- a/romancal/source_catalog/source_catalog_step.py +++ b/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 diff --git a/romancal/source_catalog/tests/test_source_catalog.py b/romancal/source_catalog/tests/test_source_catalog.py index 13580f474..60d3d9d0d 100644 --- a/romancal/source_catalog/tests/test_source_catalog.py +++ b/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