From 1b9eb0960b1f414f497ef35b6f610dae962bcfc0 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 1 Nov 2024 12:04:50 -0400 Subject: [PATCH 01/25] Add detection_image module --- romancal/multiband_catalog/detection_image.py | 149 ++++++++++++++++++ 1 file changed, 149 insertions(+) create mode 100644 romancal/multiband_catalog/detection_image.py diff --git a/romancal/multiband_catalog/detection_image.py b/romancal/multiband_catalog/detection_image.py new file mode 100644 index 000000000..03f69b08c --- /dev/null +++ b/romancal/multiband_catalog/detection_image.py @@ -0,0 +1,149 @@ +import logging +import warnings + +import numpy as np +from astropy.convolution import convolve_fft +from astropy.utils.exceptions import AstropyUserWarning +from photutils.segmentation import make_2dgaussian_kernel + +from romancal.datamodels import ModelLibrary + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + + +def make_gaussian_kernel(fwhm): + """ + Make a normalized 2D circular Gaussian kernel. + + The kernel will have odd sizes in both X and Y, be centered in the + central pixel, and normalized to sum to 1. + + Parameters + ---------- + fwhm : float + The full-width at half-maximum (FWHM) in pixels of the 2D + Gaussian kernel. + + Returns + ------- + kernel : `astropy.convolution.Kernel2D` + The output Gaussian kernel. + """ + size = np.ceil(3 * fwhm).astype(int) + size = size + 1 if size % 2 == 0 else size # make size be odd + kernel = make_2dgaussian_kernel(fwhm, size=size) # sums to 1 + return kernel.array + + +def make_det_image(library, kernel_fwhm): + """ + Make a detection image from a library of models. + + The detection image is the weighted sum of the input models, where + the weights are SED weight divided by the read noise variance. + + Parameters + ---------- + library : `romancal.datamodels.ModelLibrary` + The input library of models. + + kernel_fwhm : float + The full-width at half-maximum (FWHM) in pixels of the 2D + Gaussian kernel used to smooth the detection image. + + Returns + ------- + detection_data : 2D `numpy.ndarray` + The detection image data. + + detection_error : 2D `numpy.ndarray` + The detection image (standard deviation) error. + """ + if not isinstance(library, ModelLibrary): + raise TypeError("library input must be a ModelLibrary object") + + # TODO: extend to different kernel shapes beyond Gaussian? + # would require different/additional kernel parameters + kernel = make_gaussian_kernel(kernel_fwhm) + + log.info(f"Making detection image with kernel FWHM={kernel_fwhm}") + + with library: + detection_data = 0.0 + detection_var = 0.0 + wht_sum = 0.0 + for i, model in enumerate(library): + # TODO: SED weights to be defined in the asn file for each + # input filter image + try: + sed_weight = library.asn["products"][0]["members"][i]["sed_weight"] + except KeyError: + sed_weight = 1.0 + + log.info( + f"Processing model {model.meta.filename}: " + f"filter={model.meta.basic.optical_element}, {sed_weight=}" + ) + + # Ideally the weights should be the inverse variance of + # all sources of noise except the source Poisson noise. The + # closest approximation we have is var_rnoise (read noise + # variance) -- var_rnoise is also used as the weighting in the + # resample step. Note that model.err is the total error, which + # includes source Poisson noise. + wht = sed_weight / model.var_rnoise # inverse variance + wht_sum += wht + + coverage_mask = np.isnan(model.err) + with warnings.catch_warnings(): + # Suppress warnings about any NaNs in the data + warnings.filterwarnings(action="ignore", category=AstropyUserWarning) + + detection_data += convolve_fft( + wht * model.data, kernel, mask=coverage_mask + ) + detection_var += convolve_fft( + wht**2 * model.var_rnoise, kernel, mask=coverage_mask + ) + + library.shelve(model, modify=False) + + detection_data /= wht_sum + detection_error = np.sqrt(detection_var) / wht_sum # std dev error + + return detection_data, detection_error + + +def make_detection_image(library, kernel_fwhms): + """ + Make a detection image from a library of models. + + The detection image is maximum of the detection images from each of + the kernels. + + Parameters + ---------- + library : `romancal.datamodels.ModelLibrary` + The input library of models. + + kernel_fwhms : list of float + The full-width at half-maximum (FWHM) in pixels of the 2D + Gaussian kernels used to smooth the detection image. + + Returns + ------- + detection_data : 2D `numpy.ndarray` + The detection image data. + + detection_error : 2D `numpy.ndarray` + The detection image (standard deviation) error. + """ + log.info("Making detection image") + det_img = 0 + det_err = 0 + for kernel_fwhm in kernel_fwhms: + img, err = make_det_image(library, kernel_fwhm) + det_img = np.maximum(det_img, img) + det_err = np.maximum(det_err, err) + return det_img, det_err From 581dfd27ff193e68f618e5d456a89484d03d9e5c Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 1 Nov 2024 12:05:24 -0400 Subject: [PATCH 02/25] Add catalog utils --- romancal/multiband_catalog/utils.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 romancal/multiband_catalog/utils.py diff --git a/romancal/multiband_catalog/utils.py b/romancal/multiband_catalog/utils.py new file mode 100644 index 000000000..8b2cf4f97 --- /dev/null +++ b/romancal/multiband_catalog/utils.py @@ -0,0 +1,19 @@ +from astropy.table import Table + + +def update_colnames(table, prefix): + if not isinstance(table, Table): + raise ValueError("table must be an astropy.table.Table object") + if not isinstance(prefix, str): + raise ValueError("prefix must be a string") + + for col in table.colnames: + if ( + col.startswith("aper") + or col.startswith("CI_") + or col.startswith("isophotal_flux") + or col.startswith("kron_flux") + or "_psf" in col + ): + table.rename_column(col, prefix + col) + return table From 07f7eb3ae5195c23a3c7b62bd6a98dc58d3941f3 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 1 Nov 2024 12:07:47 -0400 Subject: [PATCH 03/25] Clean up source_catalog_step --- romancal/source_catalog/source_catalog_step.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/romancal/source_catalog/source_catalog_step.py b/romancal/source_catalog/source_catalog_step.py index e1a230b2f..1446a4e6e 100644 --- a/romancal/source_catalog/source_catalog_step.py +++ b/romancal/source_catalog/source_catalog_step.py @@ -122,7 +122,10 @@ def process(self, step_input): if segment_img is None: # no sources found source_catalog_model.source_catalog = Table() else: - ci_star_thresholds = (self.ci1_star_threshold, self.ci2_star_threshold) + ci_star_thresholds = ( + self.ci1_star_threshold, + self.ci2_star_threshold, + ) catobj = RomanSourceCatalog( model, segment_img, @@ -190,7 +193,10 @@ def save_base_results(self, segment_img, source_catalog_model): # save the source catalog self.save_model( - source_catalog_model, output_file=output_filename, suffix="cat", force=True + source_catalog_model, + output_file=output_filename, + suffix="cat", + force=True, ) From 4ec675b524fdeaccb1e46470512ac14db36803ed Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 1 Nov 2024 12:07:56 -0400 Subject: [PATCH 04/25] Update RomanSourceCatalog for multiband catalog --- romancal/source_catalog/source_catalog.py | 83 +++++++++++++++++++---- 1 file changed, 69 insertions(+), 14 deletions(-) diff --git a/romancal/source_catalog/source_catalog.py b/romancal/source_catalog/source_catalog.py index b22a5fe07..612a6ad72 100644 --- a/romancal/source_catalog/source_catalog.py +++ b/romancal/source_catalog/source_catalog.py @@ -46,7 +46,7 @@ class RomanSourceCatalog: where sources are marked by different positive integer values. A value of zero is reserved for the background. - convolved_data : data : 2D `~numpy.ndarray` + convolved_data : 2D `~numpy.ndarray` or `None` The 2D array used to calculate the source centroid and shape measurements. @@ -70,6 +70,22 @@ class RomanSourceCatalog: This is needed to calculate the DAOFind sharpness and roundness properties (DAOFind uses a special kernel that sums to zero). + fit_psf : bool + Whether to fit a PSF model to the sources. + + detection_cat : `None` or `~photutils.segmentation.SourceCatalog`, optional + A `~photutils.segmentation.SourceCatalog` object for the + detection image. The segmentation image used to create + the detection catalog must be the same one input to + ``segment_image``. If input, then the detection catalog source + centroids and morphological/shape properties will be returned + instead of calculating them from the input ``data``. The + detection catalog centroids and shape properties will also be + used to perform aperture photometry (i.e., circular and Kron). + + flux_unit : str, optional + The unit of the flux density. Default is 'uJy'. + Notes ----- ``model.err`` is assumed to be the total error array corresponding @@ -87,6 +103,7 @@ def __init__( ci_star_thresholds, kernel_fwhm, fit_psf, + detection_cat=None, flux_unit="uJy", ): @@ -99,6 +116,7 @@ def __init__( self.aperture_params = aperture_params self.kernel_sigma = kernel_fwhm * gaussian_fwhm_to_sigma self.fit_psf = fit_psf + self.detection_cat = detection_cat if len(ci_star_thresholds) != 2: raise ValueError("ci_star_thresholds must contain only 2 items") @@ -123,6 +141,9 @@ def __init__( u.Unit(self.flux_unit) ) + # needed for Kron photometry + self.segm_sourcecat = None + @lazyproperty def _pixscale_angle(self): ysize, xsize = self.model.data.shape @@ -173,7 +194,8 @@ def convert_l2_to_sb(self): # 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 + if self.convolved_data is not None: + self.convolved_data *= self.l2_conv_factor def convert_sb_to_flux_density(self): """ @@ -191,7 +213,8 @@ def convert_sb_to_flux_density(self): self.convolved_data *= self.sb_to_flux.value self.model["data"] <<= self.sb_to_flux.unit self.model["err"] <<= self.sb_to_flux.unit - self.convolved_data <<= self.sb_to_flux.unit + if self.convolved_data is not None: + self.convolved_data <<= self.sb_to_flux.unit def convert_flux_to_abmag(self, flux, flux_err): """ @@ -239,6 +262,8 @@ def segment_colnames(self): # desc["isophotal_abmag"] = "Isophotal AB magnitude" # desc["isophotal_abmag_err"] = "Isophotal AB magnitude error" desc["isophotal_area"] = "Isophotal area" + desc["kron_flux"] = "Kron flux" + desc["kron_flux_err"] = "Kron flux error" desc["semimajor_sigma"] = ( "1-sigma standard deviation along the semimajor axis of the 2D Gaussian function that has the same second-order central moments as the source" ) @@ -278,13 +303,19 @@ def set_segment_properties(self): The results are set as dynamic attributes on the class instance. """ + if self.detection_cat is not None: + detection_cat = self.detection_cat.segm_sourcecat + else: + detection_cat = None segm_cat = SourceCatalog( self.model.data, self.segment_img, convolved_data=self.convolved_data, error=self.model.err, wcs=self.wcs, + detection_cat=detection_cat, ) + self.segm_sourcecat = segm_cat self.meta.update(segm_cat.meta) @@ -293,6 +324,7 @@ def set_segment_properties(self): prop_names["isophotal_flux"] = "segment_flux" prop_names["isophotal_flux_err"] = "segment_fluxerr" prop_names["isophotal_area"] = "area" + prop_names["kron_flux_err"] = "kron_fluxerr" # dynamically set the attributes for column in self.segment_colnames: @@ -308,10 +340,17 @@ def set_segment_properties(self): @lazyproperty def _xypos(self): """ - The (x, y) source positions, defined from the segmentation - image. + The (x, y) source positions. + + If a detection catalog is input, then the detection catalog + centroids are used. Otherwise, the segment catalog centroids are + used. """ - return np.transpose((self.xcentroid, self.ycentroid)) + if self.detection_cat is not None: + xycen = (self.detection_cat.xcentroid, self.detection_cat.ycentroid) + else: + xycen = (self.xcentroid, self.ycentroid) + return np.transpose(xycen) @lazyproperty def _xypos_aper(self): @@ -963,12 +1002,27 @@ def colnames(self): """ The column name order for the final source catalog. """ - colnames = self.segment_colnames[:4] - colnames.extend(self.aperture_colnames) - colnames.extend(self.extras_colnames) - colnames.extend(self.segment_colnames[4:]) - if self.fit_psf: - colnames.extend(self.psf_photometry_colnames) + if self.detection_cat is None: + colnames = self.segment_colnames[:4] + colnames.extend(self.aperture_colnames) + colnames.extend(self.extras_colnames) + colnames.extend(self.segment_colnames[4:]) + if self.fit_psf: + colnames.extend(self.psf_photometry_colnames) + + else: + # NOTE: label is need to join the tables + colnames = [ + "label", + "isophotal_flux", + "isophotal_flux_err", + "kron_flux", + "kron_flux_err", + ] + colnames.extend(self.aperture_colnames) + if self.fit_psf: + colnames.extend(self.psf_photometry_colnames) + return colnames def _update_metadata(self): @@ -1123,11 +1177,12 @@ def do_psf_photometry(self) -> None: ) log.info("Fitting a PSF model to sources for improved astrometric precision.") + xinit, yinit = np.transpose(self._xypos) psf_photometry_table, photometry = psf.fit_psf_to_image_model( image_model=self.model, psf_model=gridded_psf_model, - x_init=self.xcentroid, - y_init=self.ycentroid, + x_init=xinit, + y_init=yinit, exclude_out_of_bounds=True, ) From 2a1d6f4f284ac48a23cc1f1806faa33140f4a2aa Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 1 Nov 2024 12:08:09 -0400 Subject: [PATCH 05/25] Add multiband catalog step --- romancal/multiband_catalog/__init__.py | 3 + .../multiband_catalog_step.py | 210 ++++++++++++++++++ romancal/multiband_catalog/tests/__init__.py | 0 .../tests/test_multiband_catalog.py | 0 4 files changed, 213 insertions(+) create mode 100644 romancal/multiband_catalog/__init__.py create mode 100644 romancal/multiband_catalog/multiband_catalog_step.py create mode 100644 romancal/multiband_catalog/tests/__init__.py create mode 100644 romancal/multiband_catalog/tests/test_multiband_catalog.py diff --git a/romancal/multiband_catalog/__init__.py b/romancal/multiband_catalog/__init__.py new file mode 100644 index 000000000..9b14307ab --- /dev/null +++ b/romancal/multiband_catalog/__init__.py @@ -0,0 +1,3 @@ +from .multiband_catalog_step import MultibandCatalogStep + +__all__ = ["MultibandCatalogStep"] diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py new file mode 100644 index 000000000..ce29d5c18 --- /dev/null +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -0,0 +1,210 @@ +""" +Module for the multiband source catalog step. +""" + +import logging + +import numpy as np +from astropy.table import Table, join +from roman_datamodels import datamodels, maker_utils + +from romancal.datamodels import ModelLibrary +from romancal.multiband_catalog.detection_image import make_detection_image +from romancal.multiband_catalog.utils import update_colnames +from romancal.source_catalog.background import RomanBackground +from romancal.source_catalog.detection import make_segmentation_image +from romancal.source_catalog.source_catalog import RomanSourceCatalog +from romancal.stpipe import RomanStep + +__all__ = ["MultibandCatalogStep"] + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + + +class MultibandCatalogStep(RomanStep): + """ + Create a multiband catalog of sources including photometry and basic + shape measurements. + + Parameters + ----------- + input : str or `ModelLibrary` + Path to an ASDF file or a `ModelLibrary`. + """ + + class_alias = "multiband_catalog" + reference_file_types = [] + + spec = """ + bkg_boxsize = integer(default=1000) # background mesh box size in pixels + kernel_fwhms = float_list(default=None) # Gaussian kernel FWHM in pixels + snr_threshold = float(default=3.0) # per-pixel SNR threshold above the bkg + npixels = integer(default=25) # min number of pixels in source + deblend = boolean(default=False) # deblend sources? + aperture_ee1 = integer(default=30) # aperture encircled energy 1 + aperture_ee2 = integer(default=50) # aperture encircled energy 2 + aperture_ee3 = integer(default=70) # aperture encircled energy 3 + ci1_star_threshold = float(default=2.0) # CI 1 star threshold + ci2_star_threshold = float(default=1.8) # CI 2 star threshold + suffix = string(default='cat') # Default suffix for output files + fit_psf = boolean(default=True) # fit source PSFs for accurate astrometry? + """ + + def process(self, library): + if isinstance(library, str): + library = ModelLibrary(library) + if not isinstance(library, ModelLibrary): + raise TypeError("library input must be a ModelLibrary object") + + # TODO: sensible defaults + if self.kernel_fwhms is None: + self.kernel_fwhms = [2.0, 20.0] + + # NOTE: I'm assuming here that all of the input images have been + # background subtracted, have the same shape, and are pixel + # aligned. + # TODO: Do we need a separate background subtraction step + # prior to this one? + + cat_model = datamodels.MosaicSourceCatalogModel + source_catalog_model = maker_utils.mk_datamodel(cat_model) + + # define the aperture parameters for the source catalog + # based on the input encircled energy fractions + # TODO: where will these values be defined? + # source_catalog will ultimately get these filter-dependent + # values from a reference file based on EE values; + # do we want filter-dependent aperture parameters for the + # multiband catalog? + aperture_ee = (self.aperture_ee1, self.aperture_ee2, self.aperture_ee3) + aperture_params = { + "aperture_radii": np.array((1.0, 2.0, 3.0)), + "aperture_corrections": np.array((1.0, 1.0, 1.0)), + "aperture_ee": aperture_ee, + "bkg_aperture_inner_radius": 5.0, + "bkg_aperture_outer_radius": 10.0, + } + + # TODO: do we want to save the det_img and det_err? + det_img, det_err = make_detection_image(library, self.kernel_fwhms) + + # estimate background rms from detection image to calculate a + # threshold for source detection + mask = np.isnan(det_img) + coverage_mask = np.isnan(det_err) + bkg = RomanBackground( + det_img, + box_size=self.bkg_boxsize, + mask=mask, + coverage_mask=coverage_mask, + ) + bkg_rms = bkg.background_rms + + segment_img = make_segmentation_image( + det_img, + snr_threshold=self.snr_threshold, + npixels=self.npixels, + bkg_rms=bkg_rms, + deblend=self.deblend, + mask=coverage_mask, + ) + + if segment_img is None: # no sources found + source_catalog_model.source_catalog = Table() + else: + ci_star_thresholds = ( + self.ci1_star_threshold, + self.ci2_star_threshold, + ) + + # this is needed for the DAOFind sharpness and roundness + # properties; are these needed for the Roman source catalog? + star_kernel_fwhm = np.min(self.kernel_fwhms) # ?? + + det_model = maker_utils.mk_datamodel(datamodels.MosaicModel) + det_model.data = det_img + det_model.err = det_err + + # TODO: this is a temporary solution to get model attributes + # currently needed in RomanSourceCatalog + with library: + model = library.borrow(0) + det_model.weight = model.weight + det_model.meta = model.meta + library.shelve(model, modify=False) + + det_catobj = RomanSourceCatalog( + det_model, + segment_img, + det_img, + aperture_params, + ci_star_thresholds, + star_kernel_fwhm, + self.fit_psf, + detection_cat=None, + ) + det_cat = update_colnames(det_catobj.catalog, prefix="det_") + + # loop over each image + with library: + for model in enumerate(library): + catobj = RomanSourceCatalog( + model, + segment_img, + None, + aperture_params, + ci_star_thresholds, + star_kernel_fwhm, + self.fit_psf, + detection_cat=det_catobj, + ) + + filter_name = model.meta.basic.optical_element + prefix = f"{filter_name}_" + cat = update_colnames(catobj.catalog, prefix=prefix) + cat.meta = None # temporary + # outer join prevents empty table if any columns have + # the same name but different values (e.g., repeated + # filter names) + det_cat = join(det_cat, cat, keys="label", join_type="outer") + library.shelve(model, modify=False) + + # put the resulting catalog in the model + source_catalog_model.source_catalog = det_cat + + # save the segmentation image and multiband catalog + # TODO: I noticed that the catalog is saved twice; + # once here and once when the step returns + self.save_base_results(segment_img, source_catalog_model) + + return source_catalog_model + + def save_base_results(self, segment_img, source_catalog_model): + # save the segmentation map and source catalog + output_filename = ( + self.output_file + if self.output_file is not None + else source_catalog_model.meta.filename + ) + + seg_model = datamodels.MosaicSegmentationMapModel + segmentation_model = maker_utils.mk_datamodel(seg_model) + for key in segmentation_model.meta.keys(): + segmentation_model.meta[key] = source_catalog_model.meta[key] + + if segment_img is not None: + segmentation_model.data = segment_img.data.astype(np.uint32) + self.save_model( + segmentation_model, + output_file=output_filename, + 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/multiband_catalog/tests/__init__.py b/romancal/multiband_catalog/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/romancal/multiband_catalog/tests/test_multiband_catalog.py b/romancal/multiband_catalog/tests/test_multiband_catalog.py new file mode 100644 index 000000000..e69de29bb From 3a4053146d5a41260b3cc9d1fffbd637bf56c9a8 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 1 Nov 2024 12:27:58 -0400 Subject: [PATCH 06/25] Add multiband_catalog step to towncrier --- pyproject.toml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 0a64a84bb..8768f9b48 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -357,3 +357,8 @@ showcontent = true directory = "source_catalog" name = "``source_catalog``" showcontent = true + +[[tool.towncrier.type]] +directory = "multiband_catalog" +name = "``multiband_catalog``" +showcontent = true From 3eb61f8af2c3eb5778f2b4750681d000a0bf864f Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 1 Nov 2024 14:42:42 -0400 Subject: [PATCH 07/25] Add multibandcatalogstep suffix --- romancal/lib/suffix.py | 1 + 1 file changed, 1 insertion(+) diff --git a/romancal/lib/suffix.py b/romancal/lib/suffix.py index 05925d409..6df60bac4 100644 --- a/romancal/lib/suffix.py +++ b/romancal/lib/suffix.py @@ -99,6 +99,7 @@ "refpixstep", "resamplestep", "sourcecatalogstep", + "multibandcatalogstep", } From f91209c51bd661c9a3084dea256b18ac5f94405d Mon Sep 17 00:00:00 2001 From: Eddie Schlafly Date: Tue, 12 Nov 2024 09:56:55 -0500 Subject: [PATCH 08/25] updates for B16 tests. --- .../multiband_catalog/multiband_catalog_step.py | 5 +++-- romancal/regtest/test_catalog.py | 4 ++++ romancal/scripts/make_regtestdata.sh | 13 ++++++++++--- romancal/step.py | 2 ++ romancal/tests/dms_requirement_tests.json | 9 +++++++++ 5 files changed, 28 insertions(+), 5 deletions(-) diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py index ce29d5c18..1a864eaa4 100644 --- a/romancal/multiband_catalog/multiband_catalog_step.py +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -69,6 +69,7 @@ def process(self, library): cat_model = datamodels.MosaicSourceCatalogModel source_catalog_model = maker_utils.mk_datamodel(cat_model) + source_catalog_model.meta.filename = library.asn['products'][0]['name'] # define the aperture parameters for the source catalog # based on the input encircled energy fractions @@ -102,7 +103,7 @@ def process(self, library): bkg_rms = bkg.background_rms segment_img = make_segmentation_image( - det_img, + det_img - np.nanmedian(det_img), snr_threshold=self.snr_threshold, npixels=self.npixels, bkg_rms=bkg_rms, @@ -148,7 +149,7 @@ def process(self, library): # loop over each image with library: - for model in enumerate(library): + for i, model in enumerate(library): catobj = RomanSourceCatalog( model, segment_img, diff --git a/romancal/regtest/test_catalog.py b/romancal/regtest/test_catalog.py index c3c055b30..a6e630119 100644 --- a/romancal/regtest/test_catalog.py +++ b/romancal/regtest/test_catalog.py @@ -12,6 +12,7 @@ @pytest.fixture( scope="module", params=[ + "r0099101001001001001_r274dp63x31y81_prompt_F158_i2d.asdf", "r0099101001001001001_F158_visit_i2d.asdf", "r0000101001001001001_0001_wfi01_cal.asdf", ], @@ -53,6 +54,9 @@ def fields(catalog): ( "ra_centroid", # DMS374 positions on ICRF "dec_centroid", # DMS374 positions on ICRF + "aper30_flux", # DMS399 aperture fluxes + "aper50_flux", # DMS399 aperture fluxes + "aper70_flux", # DMS399 aperture fluxes "aper_total_flux", # DMS375 fluxes "is_extended", # DMS376 type of source "aper_total_flux_err", # DMS386 flux uncertainties diff --git a/romancal/scripts/make_regtestdata.sh b/romancal/scripts/make_regtestdata.sh index a6031dd35..f3e97d542 100644 --- a/romancal/scripts/make_regtestdata.sh +++ b/romancal/scripts/make_regtestdata.sh @@ -202,17 +202,24 @@ cp ${l3name}_i2d.asdf $outdir/roman-pipeline/dev/truth/WFI/image/ strun romancal.step.SourceCatalogStep ${l3name}_i2d.asdf cp ${l3name}_cat.asdf $outdir/roman-pipeline/dev/truth/WFI/image/ - # L2 catalog strun romancal.step.SourceCatalogStep r0000101001001001001_0001_wfi01_cal.asdf cp r0000101001001001001_0001_wfi01_cat.asdf $outdir/roman-pipeline/dev/truth/WFI/image/ - -l3name="r0099101001001001001_F158_visit_r274dp63x31y81" +# L3 on skycell +l3name="r0099101001001001001_r274dp63x31y81_prompt_F158" asn_from_list --product-name=$l3name r0000101001001001001_0001_wfi01_cal.asdf r0000101001001001001_0002_wfi01_cal.asdf r0000101001001001001_0003_wfi01_cal.asdf -o L3_m1_asn.json strun roman_mos L3_m1_asn.json cp ${l3name}_i2d.asdf $outdir/roman-pipeline/dev/truth/WFI/image/ +# L3 skycell catalog +strun romancal.step.SourceCatalogStep ${l3name}_i2d.asdf +cp ${l3name}_cat.asdf $outdir/roman-pipeline/dev/truth/WFI/image/ + +# multiband catalog +asn_from_list --product-name=${l3name}_mbcat ${l3name}_i2d.asdf -o L3_skycell_mbcat_asn.json +strun romancal.step.MultibandCatalogStep L3_skycell_mbcat_asn.json --deblend True + # tests passing suffix to the pipeline strun roman_elp r0000101001001001001_0001_wfi01_uncal.asdf --steps.tweakreg.skip=True --suffix=star cp r0000101001001001001_0001_wfi01_star.asdf $outdir/roman-pipeline/dev/truth/WFI/image/ diff --git a/romancal/step.py b/romancal/step.py index 8ddc7362d..b2b9a8ff2 100644 --- a/romancal/step.py +++ b/romancal/step.py @@ -19,6 +19,7 @@ from .skymatch.skymatch_step import SkyMatchStep from .source_catalog.source_catalog_step import SourceCatalogStep from .source_detection.source_detection_step import SourceDetectionStep +from .multiband_catalog.multiband_catalog_step import MultibandCatalogStep from .tweakreg.tweakreg_step import TweakRegStep __all__ = [ @@ -38,5 +39,6 @@ "SkyMatchStep", "SourceDetectionStep", "SourceCatalogStep", + "MultibandCatalogStep", "TweakRegStep", ] diff --git a/romancal/tests/dms_requirement_tests.json b/romancal/tests/dms_requirement_tests.json index 2d6588492..f3069cf82 100644 --- a/romancal/tests/dms_requirement_tests.json +++ b/romancal/tests/dms_requirement_tests.json @@ -108,6 +108,15 @@ "DMS387": [ "romancal.regtest.test_catalog.test_has_field" ], + "DMS391": [ + "romancal.regtest.test_multiband_catalog.test_multiband_catalog", + ], + "DMS393": [ + "romancal.regtest.test_multiband_catalog.test_multiband_catalog", + ], + "DMS399": [ + "romancal.regtest.test_catalog.test_has_field" + ], "DMS400": [ "romancal.regtest.test_mos_pipeline.test_steps_ran", "romancal.regtest.test_mos_pipeline.test_added_background", From f94148d0a3d039167e3f5472bd7f82befcc9de26 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Mon, 11 Nov 2024 14:15:43 -0500 Subject: [PATCH 09/25] Add background subtraction --- romancal/multiband_catalog/background.py | 34 +++++++++++++++++++ .../multiband_catalog_step.py | 20 ++++++----- 2 files changed, 46 insertions(+), 8 deletions(-) create mode 100644 romancal/multiband_catalog/background.py diff --git a/romancal/multiband_catalog/background.py b/romancal/multiband_catalog/background.py new file mode 100644 index 000000000..ff2683c7e --- /dev/null +++ b/romancal/multiband_catalog/background.py @@ -0,0 +1,34 @@ +import numpy as np +from roman_datamodels.datamodels import ImageModel, MosaicModel + +from romancal.datamodels import ModelLibrary +from romancal.source_catalog.background import RomanBackground + + +def subtract_background(model, box_size=1000): + if not isinstance(model, (ImageModel, MosaicModel)): + raise ValueError("The input model must be an ImageModel or MosaicModel.") + + # Subtract the background + mask = np.isnan(model.data) + coverage_mask = np.isnan(model.err) + bkg = RomanBackground( + model.data, + box_size=box_size, + mask=mask, + coverage_mask=coverage_mask, + ) + model.data -= bkg.background + return model + + +def subtract_background_library(library, box_size=1000): + if not isinstance(library, ModelLibrary): + raise TypeError("library input must be a ModelLibrary object") + + with library: + for model in library: + model = subtract_background(model, box_size=box_size) + library.shelve(model) + + return library diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py index 1a864eaa4..2d03834ce 100644 --- a/romancal/multiband_catalog/multiband_catalog_step.py +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -9,6 +9,7 @@ from roman_datamodels import datamodels, maker_utils from romancal.datamodels import ModelLibrary +from romancal.multiband_catalog.background import subtract_background_library from romancal.multiband_catalog.detection_image import make_detection_image from romancal.multiband_catalog.utils import update_colnames from romancal.source_catalog.background import RomanBackground @@ -29,8 +30,10 @@ class MultibandCatalogStep(RomanStep): Parameters ----------- - input : str or `ModelLibrary` - Path to an ASDF file or a `ModelLibrary`. + input : str or `~romancal.datamodels.ModelLibrary` + Path to an ASDF file or a `~romancal.datamodels.ModelLibrary` + that contains `~roman_datamodels.datamodels.MosaicImageModel` + models. """ class_alias = "multiband_catalog" @@ -52,6 +55,8 @@ class MultibandCatalogStep(RomanStep): """ def process(self, library): + # All input MosaicImages in the ModelLibrary are assumed to have + # the same shape and be pixel aligned. if isinstance(library, str): library = ModelLibrary(library) if not isinstance(library, ModelLibrary): @@ -61,11 +66,8 @@ def process(self, library): if self.kernel_fwhms is None: self.kernel_fwhms = [2.0, 20.0] - # NOTE: I'm assuming here that all of the input images have been - # background subtracted, have the same shape, and are pixel - # aligned. - # TODO: Do we need a separate background subtraction step - # prior to this one? + # NOTE: I'm assuming here that all of the input images have + # have the same shape and are pixel aligned. cat_model = datamodels.MosaicSourceCatalogModel source_catalog_model = maker_utils.mk_datamodel(cat_model) @@ -87,6 +89,8 @@ def process(self, library): "bkg_aperture_outer_radius": 10.0, } + library = subtract_background_library(library, self.bkg_boxsize) + # TODO: do we want to save the det_img and det_err? det_img, det_err = make_detection_image(library, self.kernel_fwhms) @@ -149,7 +153,7 @@ def process(self, library): # loop over each image with library: - for i, model in enumerate(library): + for model in library: catobj = RomanSourceCatalog( model, segment_img, From d898f49cadf7cb8d8f3e05ba1c7a15f955c4cad6 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Mon, 11 Nov 2024 23:22:24 -0500 Subject: [PATCH 10/25] Remove columns with det_ prefix --- romancal/multiband_catalog/utils.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/romancal/multiband_catalog/utils.py b/romancal/multiband_catalog/utils.py index 8b2cf4f97..39704ff63 100644 --- a/romancal/multiband_catalog/utils.py +++ b/romancal/multiband_catalog/utils.py @@ -2,6 +2,26 @@ def update_colnames(table, prefix): + """ + Update column names in an astropy table by adding a prefix to + photometry-related columns. + + If the prefix is "det_", the function will remove these columns. + + Parameters + ---------- + table : astropy.table.Table + The table to update. + + prefix : str + The prefix to add to the photometry-related column names. If + "det_", the columns will be removed. + + Returns + ------- + astropy.table.Table + The updated table. + """ if not isinstance(table, Table): raise ValueError("table must be an astropy.table.Table object") if not isinstance(prefix, str): @@ -15,5 +35,9 @@ def update_colnames(table, prefix): or col.startswith("kron_flux") or "_psf" in col ): - table.rename_column(col, prefix + col) + if prefix == "det_": + table.remove_column(col) + else: + table.rename_column(col, prefix + col) + return table From 5e07605adc1b8dcc935a0232ed970ff71d05ccb4 Mon Sep 17 00:00:00 2001 From: Eddie Schlafly Date: Tue, 12 Nov 2024 11:40:57 -0500 Subject: [PATCH 11/25] Actually add new regtests. --- romancal/regtest/test_multiband_catalog.py | 60 ++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 romancal/regtest/test_multiband_catalog.py diff --git a/romancal/regtest/test_multiband_catalog.py b/romancal/regtest/test_multiband_catalog.py new file mode 100644 index 000000000..6114ae32c --- /dev/null +++ b/romancal/regtest/test_multiband_catalog.py @@ -0,0 +1,60 @@ +""" Roman tests for source catalog creation """ + +import asdf +import pytest + +from romancal.stpipe import RomanStep +from romancal.multiband_catalog.multiband_catalog_step import MultibandCatalogStep + +# mark all tests in this module +pytestmark = [pytest.mark.bigdata, pytest.mark.soctests] + +fieldlist = [ + "ra_centroid", # DMS374 positions on ICRF + "dec_centroid", # DMS374 positions on ICRF + "F158_aper30_flux", # DMS399 aperture fluxes + "F158_aper50_flux", # DMS399 aperture fluxes + "F158_aper70_flux", # DMS399 aperture fluxes + "F158_aper_total_flux", # DMS375 fluxes + "F158_aper_total_flux_err", # DMS386 flux uncertainties + "flags", # DMS387 dq_flags +] + + +def test_multiband_catalog(rtdata_module): + rtdata = rtdata_module + inputasnfn = "L3_skycell_mbcat_asn.json" + # note that this input association currently only has a single + # filter in it, so this is more of an existence proof for the multiband + # catalogs than a detailed test. Using only a single catalog lets us + # rely on the existing regtest files. + outputfn = "r0099101001001001001_r274dp63x31y81_prompt_cat.asdf" + rtdata.get_asn(f"WFI/image/{inputasnfn}") + rtdata.output = outputfn + rtdata.input = inputasnfn + rtdata.get_truth(f"truth/WFI/image/{outputfn}") + aftruth = asdf.open(f"truth/{outputfn}") + args = ['romancal.step.MultibandCatalogStep', inputasnfn, + '--deblend', 'True', # use deblending, DMS 393 + ] + RomanStep.from_cmdline(args) + afcat = asdf.open(outputfn) + for field in fieldlist: + assert field in afcat['roman']['source_catalog'].dtype.names + + # DMS 393: multiband catalog uses both PSF-like and extend-source-like + # kernels + assert ( + set(aftruth['roman']['source_catalog'].dtype.names) == + set(afcat['roman']['source_catalog'].dtype.names)) + # weak assertion that our truth file must at least have the same + # catalog fields as the file produced here. Exactly matching rows + # would require a lot of okifying things that aren't obviously + # the same, but we can easily check that the columns match up + # by name. + + step = MultibandCatalogStep() + step.log.info('DMS391: successfully used multiple kernels to detect ' + 'sources.') + step.log.info('DMS393: successfully used deblending to separate blended ' + 'sources.') From c6dd7c01358bc645f45b9b6909e1a7f2c740dbb6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 12 Nov 2024 15:44:26 +0000 Subject: [PATCH 12/25] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- romancal/multiband_catalog/multiband_catalog_step.py | 2 +- romancal/step.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py index 2d03834ce..aa01e1278 100644 --- a/romancal/multiband_catalog/multiband_catalog_step.py +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -71,7 +71,7 @@ def process(self, library): cat_model = datamodels.MosaicSourceCatalogModel source_catalog_model = maker_utils.mk_datamodel(cat_model) - source_catalog_model.meta.filename = library.asn['products'][0]['name'] + source_catalog_model.meta.filename = library.asn["products"][0]["name"] # define the aperture parameters for the source catalog # based on the input encircled energy fractions diff --git a/romancal/step.py b/romancal/step.py index b2b9a8ff2..e365efb4d 100644 --- a/romancal/step.py +++ b/romancal/step.py @@ -10,6 +10,7 @@ from .flux import FluxStep from .jump.jump_step import JumpStep from .linearity.linearity_step import LinearityStep +from .multiband_catalog.multiband_catalog_step import MultibandCatalogStep from .outlier_detection.outlier_detection_step import OutlierDetectionStep from .photom.photom_step import PhotomStep from .ramp_fitting.ramp_fit_step import RampFitStep @@ -19,7 +20,6 @@ from .skymatch.skymatch_step import SkyMatchStep from .source_catalog.source_catalog_step import SourceCatalogStep from .source_detection.source_detection_step import SourceDetectionStep -from .multiband_catalog.multiband_catalog_step import MultibandCatalogStep from .tweakreg.tweakreg_step import TweakRegStep __all__ = [ From ef70c97d5bf9bba6b1d89a43a6b382606624e003 Mon Sep 17 00:00:00 2001 From: Eddie Schlafly Date: Tue, 12 Nov 2024 11:43:24 -0500 Subject: [PATCH 13/25] Increase number of ids. --- romancal/regtest/test_catalog.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/romancal/regtest/test_catalog.py b/romancal/regtest/test_catalog.py index a6e630119..8f525fd4e 100644 --- a/romancal/regtest/test_catalog.py +++ b/romancal/regtest/test_catalog.py @@ -16,7 +16,7 @@ "r0099101001001001001_F158_visit_i2d.asdf", "r0000101001001001001_0001_wfi01_cal.asdf", ], - ids=["L3", "L2"], + ids=["L3", "L2", "L3skycell"], ) def run_source_catalog(rtdata_module, request): rtdata = rtdata_module From aa9e228d5e77ecc089f907f4cd7087ec29b45a36 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 12 Nov 2024 16:42:18 +0000 Subject: [PATCH 14/25] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- romancal/regtest/test_multiband_catalog.py | 29 ++++++++++++---------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/romancal/regtest/test_multiband_catalog.py b/romancal/regtest/test_multiband_catalog.py index 6114ae32c..f9c2a8694 100644 --- a/romancal/regtest/test_multiband_catalog.py +++ b/romancal/regtest/test_multiband_catalog.py @@ -3,8 +3,8 @@ import asdf import pytest -from romancal.stpipe import RomanStep from romancal.multiband_catalog.multiband_catalog_step import MultibandCatalogStep +from romancal.stpipe import RomanStep # mark all tests in this module pytestmark = [pytest.mark.bigdata, pytest.mark.soctests] @@ -34,19 +34,22 @@ def test_multiband_catalog(rtdata_module): rtdata.input = inputasnfn rtdata.get_truth(f"truth/WFI/image/{outputfn}") aftruth = asdf.open(f"truth/{outputfn}") - args = ['romancal.step.MultibandCatalogStep', inputasnfn, - '--deblend', 'True', # use deblending, DMS 393 - ] + args = [ + "romancal.step.MultibandCatalogStep", + inputasnfn, + "--deblend", + "True", # use deblending, DMS 393 + ] RomanStep.from_cmdline(args) afcat = asdf.open(outputfn) for field in fieldlist: - assert field in afcat['roman']['source_catalog'].dtype.names - + assert field in afcat["roman"]["source_catalog"].dtype.names + # DMS 393: multiband catalog uses both PSF-like and extend-source-like # kernels - assert ( - set(aftruth['roman']['source_catalog'].dtype.names) == - set(afcat['roman']['source_catalog'].dtype.names)) + assert set(aftruth["roman"]["source_catalog"].dtype.names) == set( + afcat["roman"]["source_catalog"].dtype.names + ) # weak assertion that our truth file must at least have the same # catalog fields as the file produced here. Exactly matching rows # would require a lot of okifying things that aren't obviously @@ -54,7 +57,7 @@ def test_multiband_catalog(rtdata_module): # by name. step = MultibandCatalogStep() - step.log.info('DMS391: successfully used multiple kernels to detect ' - 'sources.') - step.log.info('DMS393: successfully used deblending to separate blended ' - 'sources.') + step.log.info("DMS391: successfully used multiple kernels to detect " "sources.") + step.log.info( + "DMS393: successfully used deblending to separate blended " "sources." + ) From 4e3bcac5bdb0bb0c52a9db9a2c1f2219fbf376c0 Mon Sep 17 00:00:00 2001 From: Eddie Schlafly Date: Tue, 12 Nov 2024 13:24:59 -0500 Subject: [PATCH 15/25] Update log messages and add multiband catalog to integration.py --- romancal/regtest/test_multiband_catalog.py | 6 ++++-- romancal/stpipe/integration.py | 1 + romancal/tests/dms_requirement_tests.json | 6 +++--- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/romancal/regtest/test_multiband_catalog.py b/romancal/regtest/test_multiband_catalog.py index f9c2a8694..c637b7da2 100644 --- a/romancal/regtest/test_multiband_catalog.py +++ b/romancal/regtest/test_multiband_catalog.py @@ -57,7 +57,9 @@ def test_multiband_catalog(rtdata_module): # by name. step = MultibandCatalogStep() - step.log.info("DMS391: successfully used multiple kernels to detect " "sources.") + step.log.info("DMS391: successfully used multiple kernels to detect sources.") step.log.info( - "DMS393: successfully used deblending to separate blended " "sources." + "DMS393: successfully used deblending to separate blended sources." ) + step.log.info("DMS399: successfully tested that catalogs contain aperture " + "fluxes and uncertainties.") diff --git a/romancal/stpipe/integration.py b/romancal/stpipe/integration.py index 31094ce9e..362464c81 100644 --- a/romancal/stpipe/integration.py +++ b/romancal/stpipe/integration.py @@ -40,4 +40,5 @@ def get_steps(): ("romancal.step.TweakRegStep", "tweakreg", False), ("romancal.step.ResampleStep", "resample", False), ("romancal.step.SourceCatalogStep", "source_catalog", False), + ("romancal.step.MultibandCatalogStep", "multiband_catalog", False), ] diff --git a/romancal/tests/dms_requirement_tests.json b/romancal/tests/dms_requirement_tests.json index f3069cf82..2ac3043c5 100644 --- a/romancal/tests/dms_requirement_tests.json +++ b/romancal/tests/dms_requirement_tests.json @@ -109,13 +109,13 @@ "romancal.regtest.test_catalog.test_has_field" ], "DMS391": [ - "romancal.regtest.test_multiband_catalog.test_multiband_catalog", + "romancal.regtest.test_multiband_catalog.test_multiband_catalog" ], "DMS393": [ - "romancal.regtest.test_multiband_catalog.test_multiband_catalog", + "romancal.regtest.test_multiband_catalog.test_multiband_catalog" ], "DMS399": [ - "romancal.regtest.test_catalog.test_has_field" + "romancal.regtest.test_multiband_catalog.test_multiband_catalog" ], "DMS400": [ "romancal.regtest.test_mos_pipeline.test_steps_ran", From ebcc88e8d30925694c9ae762d02696a1fe202153 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 12 Nov 2024 18:25:19 +0000 Subject: [PATCH 16/25] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- romancal/regtest/test_multiband_catalog.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/romancal/regtest/test_multiband_catalog.py b/romancal/regtest/test_multiband_catalog.py index c637b7da2..c408c820d 100644 --- a/romancal/regtest/test_multiband_catalog.py +++ b/romancal/regtest/test_multiband_catalog.py @@ -58,8 +58,8 @@ def test_multiband_catalog(rtdata_module): step = MultibandCatalogStep() step.log.info("DMS391: successfully used multiple kernels to detect sources.") + step.log.info("DMS393: successfully used deblending to separate blended sources.") step.log.info( - "DMS393: successfully used deblending to separate blended sources." + "DMS399: successfully tested that catalogs contain aperture " + "fluxes and uncertainties." ) - step.log.info("DMS399: successfully tested that catalogs contain aperture " - "fluxes and uncertainties.") From 4066ad2d42bcc91fbbfa2a195dae4d86a6db4015 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Tue, 12 Nov 2024 12:42:56 -0500 Subject: [PATCH 17/25] Improve column renaming/removal --- .../multiband_catalog_step.py | 6 +- romancal/multiband_catalog/utils.py | 77 ++++++++++++++----- 2 files changed, 62 insertions(+), 21 deletions(-) diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py index aa01e1278..271577474 100644 --- a/romancal/multiband_catalog/multiband_catalog_step.py +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -11,7 +11,7 @@ from romancal.datamodels import ModelLibrary from romancal.multiband_catalog.background import subtract_background_library from romancal.multiband_catalog.detection_image import make_detection_image -from romancal.multiband_catalog.utils import update_colnames +from romancal.multiband_catalog.utils import prefix_colnames, remove_columns from romancal.source_catalog.background import RomanBackground from romancal.source_catalog.detection import make_segmentation_image from romancal.source_catalog.source_catalog import RomanSourceCatalog @@ -149,7 +149,7 @@ def process(self, library): self.fit_psf, detection_cat=None, ) - det_cat = update_colnames(det_catobj.catalog, prefix="det_") + det_cat = remove_columns(det_catobj.catalog) # loop over each image with library: @@ -167,7 +167,7 @@ def process(self, library): filter_name = model.meta.basic.optical_element prefix = f"{filter_name}_" - cat = update_colnames(catobj.catalog, prefix=prefix) + cat = prefix_colnames(catobj.catalog, prefix=prefix) cat.meta = None # temporary # outer join prevents empty table if any columns have # the same name but different values (e.g., repeated diff --git a/romancal/multiband_catalog/utils.py b/romancal/multiband_catalog/utils.py index 39704ff63..b9e6f5320 100644 --- a/romancal/multiband_catalog/utils.py +++ b/romancal/multiband_catalog/utils.py @@ -1,32 +1,24 @@ from astropy.table import Table -def update_colnames(table, prefix): +def get_photometry_columns(table): """ - Update column names in an astropy table by adding a prefix to - photometry-related columns. - - If the prefix is "det_", the function will remove these columns. + Get the photometry-related column names in an astropy table. Parameters ---------- table : astropy.table.Table - The table to update. - - prefix : str - The prefix to add to the photometry-related column names. If - "det_", the columns will be removed. + The table to search for photometry-related columns. Returns ------- - astropy.table.Table - The updated table. + result: list + The list of photometry-related column names. """ if not isinstance(table, Table): raise ValueError("table must be an astropy.table.Table object") - if not isinstance(prefix, str): - raise ValueError("prefix must be a string") + phot_cols = [] for col in table.colnames: if ( col.startswith("aper") @@ -35,9 +27,58 @@ def update_colnames(table, prefix): or col.startswith("kron_flux") or "_psf" in col ): - if prefix == "det_": - table.remove_column(col) - else: - table.rename_column(col, prefix + col) + phot_cols.append(col) + + return phot_cols + + +def prefix_colnames(table, prefix): + """ + Update column names in an astropy table in-place by adding a prefix + to photometry-related columns. + + Parameters + ---------- + table : `~astropy.table.Table` + The table to update. + + prefix : str + The prefix to add to the photometry-related column names. If + "det_", the columns will be removed. + + Returns + ------- + result : `~astropy.table.Table` + The updated table. + """ + if not isinstance(prefix, str): + raise ValueError("prefix must be a string") + + phot_cols = get_photometry_columns(table) + for col in table.colnames: + if col in phot_cols: + table.rename_column(col, prefix + col) + + return table + + +def remove_columns(table): + """ + Remove photometry-related columns from an astropy table in-place. + + Parameters + ---------- + table : `~astropy.table.Table` + The table to update. + + Returns + ------- + result : `~astropy.table.Table` + The updated table. + """ + phot_cols = get_photometry_columns(table) + for col in table.colnames: + if col in phot_cols: + table.remove_column(col) return table From 0e709a9ef299193bd6ce75212f93f97350df3686 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Tue, 12 Nov 2024 12:51:30 -0500 Subject: [PATCH 18/25] Include err==0 in the background coverage mask --- romancal/multiband_catalog/background.py | 34 +++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/romancal/multiband_catalog/background.py b/romancal/multiband_catalog/background.py index ff2683c7e..2f602cd37 100644 --- a/romancal/multiband_catalog/background.py +++ b/romancal/multiband_catalog/background.py @@ -6,12 +6,28 @@ def subtract_background(model, box_size=1000): + """ + Subtract the background from the input model. + + Parameters + ---------- + model : ImageModel or MosaicModel + The input model to subtract the background from. + + box_size : int + The size of the box to use for the background estimation. + + Returns + ------- + model : ImageModel or MosaicModel + The input model with the background subtracted. + """ if not isinstance(model, (ImageModel, MosaicModel)): raise ValueError("The input model must be an ImageModel or MosaicModel.") # Subtract the background mask = np.isnan(model.data) - coverage_mask = np.isnan(model.err) + coverage_mask = np.logical_or(np.isnan(model.err), model.err == 0) bkg = RomanBackground( model.data, box_size=box_size, @@ -23,6 +39,22 @@ def subtract_background(model, box_size=1000): def subtract_background_library(library, box_size=1000): + """ + Subtract the background from all models in the input library. + + Parameters + ---------- + library : ModelLibrary + The input library to subtract the background from. + + box_size : int + The size of the box to use for the background estimation. + + Returns + ------- + library : ModelLibrary + The input library with the background subtracted. + """ if not isinstance(library, ModelLibrary): raise TypeError("library input must be a ModelLibrary object") From 6511eff02d5b88bd85bb520b94319bb0e50bc74d Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Tue, 12 Nov 2024 13:32:10 -0500 Subject: [PATCH 19/25] Explicitly save only the segmentation image --- .../multiband_catalog_step.py | 50 ++++++++----------- 1 file changed, 21 insertions(+), 29 deletions(-) diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py index 271577474..cd2dfb3c7 100644 --- a/romancal/multiband_catalog/multiband_catalog_step.py +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -62,24 +62,18 @@ def process(self, library): if not isinstance(library, ModelLibrary): raise TypeError("library input must be a ModelLibrary object") - # TODO: sensible defaults - if self.kernel_fwhms is None: - self.kernel_fwhms = [2.0, 20.0] - - # NOTE: I'm assuming here that all of the input images have - # have the same shape and are pixel aligned. - cat_model = datamodels.MosaicSourceCatalogModel source_catalog_model = maker_utils.mk_datamodel(cat_model) source_catalog_model.meta.filename = library.asn["products"][0]["name"] + # TODO: sensible defaults + # TODO: redefine in terms of intrinsic FWHM + if self.kernel_fwhms is None: + self.kernel_fwhms = [2.0, 20.0] + # define the aperture parameters for the source catalog # based on the input encircled energy fractions - # TODO: where will these values be defined? - # source_catalog will ultimately get these filter-dependent - # values from a reference file based on EE values; - # do we want filter-dependent aperture parameters for the - # multiband catalog? + # TODO: define these values in RomanSourceCatalog aperture_ee = (self.aperture_ee1, self.aperture_ee2, self.aperture_ee3) aperture_params = { "aperture_radii": np.array((1.0, 2.0, 3.0)), @@ -91,7 +85,7 @@ def process(self, library): library = subtract_background_library(library, self.bkg_boxsize) - # TODO: do we want to save the det_img and det_err? + # TODO: where to save the det_img and det_err? det_img, det_err = make_detection_image(library, self.kernel_fwhms) # estimate background rms from detection image to calculate a @@ -107,7 +101,7 @@ def process(self, library): bkg_rms = bkg.background_rms segment_img = make_segmentation_image( - det_img - np.nanmedian(det_img), + det_img, snr_threshold=self.snr_threshold, npixels=self.npixels, bkg_rms=bkg_rms, @@ -123,8 +117,10 @@ def process(self, library): self.ci2_star_threshold, ) - # this is needed for the DAOFind sharpness and roundness - # properties; are these needed for the Roman source catalog? + # this is needed for the DAOStarFinder sharpness and roundness + # properties + # TODO: measure on a secondary detection image with minimal + # smoothing; same for basic shape parameters star_kernel_fwhm = np.min(self.kernel_fwhms) # ?? det_model = maker_utils.mk_datamodel(datamodels.MosaicModel) @@ -168,6 +164,8 @@ def process(self, library): filter_name = model.meta.basic.optical_element prefix = f"{filter_name}_" cat = prefix_colnames(catobj.catalog, prefix=prefix) + # TODO: what metadata do we want to keep, if any, + # from the filter catalogs? cat.meta = None # temporary # outer join prevents empty table if any columns have # the same name but different values (e.g., repeated @@ -178,15 +176,16 @@ def process(self, library): # put the resulting catalog in the model source_catalog_model.source_catalog = det_cat - # save the segmentation image and multiband catalog - # TODO: I noticed that the catalog is saved twice; - # once here and once when the step returns - self.save_base_results(segment_img, source_catalog_model) + # explicitly save the segmentation image; + # the multiband catalog is automatically saved by stpipe + self.save_segm(segment_img, source_catalog_model) return source_catalog_model - def save_base_results(self, segment_img, source_catalog_model): - # save the segmentation map and source catalog + def save_segm(self, segment_img, source_catalog_model): + """ + Save the segmentation image. + """ output_filename = ( self.output_file if self.output_file is not None @@ -206,10 +205,3 @@ 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, - ) From e48c16acec820615b24eecff633bf1627d1cbd81 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Tue, 12 Nov 2024 13:42:58 -0500 Subject: [PATCH 20/25] Temporarily change bkg_boxsize default --- romancal/multiband_catalog/multiband_catalog_step.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py index cd2dfb3c7..ea9809f6d 100644 --- a/romancal/multiband_catalog/multiband_catalog_step.py +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -40,7 +40,7 @@ class MultibandCatalogStep(RomanStep): reference_file_types = [] spec = """ - bkg_boxsize = integer(default=1000) # background mesh box size in pixels + bkg_boxsize = integer(default=100) # background mesh box size in pixels kernel_fwhms = float_list(default=None) # Gaussian kernel FWHM in pixels snr_threshold = float(default=3.0) # per-pixel SNR threshold above the bkg npixels = integer(default=25) # min number of pixels in source From e128953dcc79dc2977d6723bf4d09b0b83c6b726 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Tue, 12 Nov 2024 13:54:52 -0500 Subject: [PATCH 21/25] Convolve variance with kernel**2 --- romancal/multiband_catalog/detection_image.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/romancal/multiband_catalog/detection_image.py b/romancal/multiband_catalog/detection_image.py index 03f69b08c..bebc03ef9 100644 --- a/romancal/multiband_catalog/detection_image.py +++ b/romancal/multiband_catalog/detection_image.py @@ -104,7 +104,11 @@ def make_det_image(library, kernel_fwhm): wht * model.data, kernel, mask=coverage_mask ) detection_var += convolve_fft( - wht**2 * model.var_rnoise, kernel, mask=coverage_mask + wht**2 * model.var_rnoise, + kernel**2, + mask=coverage_mask, + normalize_kernel=False, + nan_treatment="fill", ) library.shelve(model, modify=False) From 52df7422932bb0e09a5c68c076c30ea94d60b994 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Tue, 12 Nov 2024 14:35:42 -0500 Subject: [PATCH 22/25] Add MultibandCatalogStep tests --- .../multiband_catalog_step.py | 5 +- .../tests/test_multiband_catalog.py | 130 ++++++++++++++++++ 2 files changed, 134 insertions(+), 1 deletion(-) diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py index ea9809f6d..2d0079779 100644 --- a/romancal/multiband_catalog/multiband_catalog_step.py +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -64,7 +64,10 @@ def process(self, library): cat_model = datamodels.MosaicSourceCatalogModel source_catalog_model = maker_utils.mk_datamodel(cat_model) - source_catalog_model.meta.filename = library.asn["products"][0]["name"] + try: + source_catalog_model.meta.filename = library.asn["products"][0]["name"] + except (AttributeError, KeyError): + source_catalog_model.meta.filename = "multiband_catalog" # TODO: sensible defaults # TODO: redefine in terms of intrinsic FWHM diff --git a/romancal/multiband_catalog/tests/test_multiband_catalog.py b/romancal/multiband_catalog/tests/test_multiband_catalog.py index e69de29bb..c73cd9262 100644 --- a/romancal/multiband_catalog/tests/test_multiband_catalog.py +++ b/romancal/multiband_catalog/tests/test_multiband_catalog.py @@ -0,0 +1,130 @@ +import os + +import astropy.units as u +import numpy as np +import pytest +from astropy.modeling.models import Gaussian2D +from astropy.table import Table +from roman_datamodels.datamodels import MosaicModel +from roman_datamodels.maker_utils import mk_level3_mosaic + +from romancal.datamodels import ModelLibrary +from romancal.multiband_catalog import MultibandCatalogStep + + +def make_test_image(): + g1 = Gaussian2D(121.0, 11.1, 12.2, 1.5, 1.5) + g2 = Gaussian2D(70, 65, 18, 9.2, 4.5) + g3 = Gaussian2D(111.0, 41, 42.7, 8.0, 3.0, theta=30 * u.deg) + g4 = Gaussian2D(81.0, 17, 52.7, 4, 2, theta=102 * u.deg) + g5 = Gaussian2D(107.0, 65, 71, 12, 2, theta=142 * u.deg) + g6 = Gaussian2D(50, 20, 80, 2.1, 2.1) + g7 = Gaussian2D(97.0, 85, 87.3, 4, 2, theta=-30 * u.deg) + + yy, xx = np.mgrid[0:101, 0:101] + data = ( + g1(xx, yy) + + g2(xx, yy) + + g3(xx, yy) + + g4(xx, yy) + + g5(xx, yy) + + g6(xx, yy) + + g7(xx, yy) + ).value.astype("float32") + + rng = np.random.default_rng(seed=123) + noise = rng.normal(0, 2.5, size=data.shape) + data += noise + err = data / 10.0 + + return data, err + + +@pytest.fixture +def mosaic_model(): + wfi_mosaic = mk_level3_mosaic(shape=(101, 101)) + model = MosaicModel(wfi_mosaic) + data, err = make_test_image() + model.data = data + model.err = err + model.var_rnoise = err**2 + model.weight = 1.0 / err + return model + + +@pytest.fixture +def library_model(mosaic_model): + model2 = mosaic_model.copy() + model2.meta.basic.optical_element = "F184" + return ModelLibrary([mosaic_model, model2]) + + +@pytest.mark.webbpsf +@pytest.mark.parametrize( + "snr_threshold, npixels, save_results", + ( + (3, 10, True), + (10, 10, False), + ), +) +def test_multiband_catalog( + library_model, snr_threshold, npixels, save_results, tmp_path +): + os.chdir(tmp_path) + step = MultibandCatalogStep() + + result = step.call( + library_model, + bkg_boxsize=50, + snr_threshold=snr_threshold, + npixels=npixels, + fit_psf=False, + save_results=save_results, + ) + + cat = result.source_catalog + + assert isinstance(cat, Table) + + columns = [ + "label", + "xcentroid", + "ycentroid", + "ra_centroid", + "dec_centroid", + "is_extended", + "sharpness", + "roundness", + "nn_label", + "nn_dist", + "isophotal_area", + "semimajor_sigma", + "semiminor_sigma", + "ellipticity", + "orientation", + "sky_orientation", + "ra_bbox_ll", + "dec_bbox_ll", + "ra_bbox_ul", + "dec_bbox_ul", + "ra_bbox_lr", + "dec_bbox_lr", + "ra_bbox_ur", + "dec_bbox_ur", + ] + + assert len(cat) == 7 + + if len(cat) > 0: + for col in columns: + assert col in cat.colnames + assert np.min(cat["xcentroid"]) > 0.0 + assert np.min(cat["ycentroid"]) > 0.0 + assert np.max(cat["xcentroid"]) < 100.0 + assert np.max(cat["ycentroid"]) < 100.0 + + for phottype in ("isophotal", "kron", "aper30", "aper_total"): + for filt in ("F158", "F184"): + colname = f"{filt}_{phottype}_flux" + assert colname in cat.colnames + assert colname + "_err" in cat.colnames From e5ff76506924cbc3a4bf32c2670612feb7c26e21 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 1 Nov 2024 12:21:18 -0400 Subject: [PATCH 23/25] Add changelog entry --- changes/1485.multiband_catalog.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 changes/1485.multiband_catalog.rst diff --git a/changes/1485.multiband_catalog.rst b/changes/1485.multiband_catalog.rst new file mode 100644 index 000000000..828e990b7 --- /dev/null +++ b/changes/1485.multiband_catalog.rst @@ -0,0 +1 @@ +Added a pipeline step to create a multiband catalog from L3 images. From ed04e7c9daba0b7dd5599a46139164d19d011f60 Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Thu, 14 Nov 2024 14:22:22 -0500 Subject: [PATCH 24/25] Check if convolved_data is None --- romancal/source_catalog/source_catalog.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/romancal/source_catalog/source_catalog.py b/romancal/source_catalog/source_catalog.py index 612a6ad72..f76c920c7 100644 --- a/romancal/source_catalog/source_catalog.py +++ b/romancal/source_catalog/source_catalog.py @@ -209,11 +209,11 @@ def convert_sb_to_flux_density(self): # 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.sb_to_flux.value - self.model["err"] *= self.sb_to_flux.value - self.convolved_data *= self.sb_to_flux.value self.model["data"] <<= self.sb_to_flux.unit + self.model["err"] *= self.sb_to_flux.value self.model["err"] <<= self.sb_to_flux.unit if self.convolved_data is not None: + self.convolved_data *= self.sb_to_flux.value self.convolved_data <<= self.sb_to_flux.unit def convert_flux_to_abmag(self, flux, flux_err): From cc0202eaf1c1bc32d2ee0823c8f1abd94aecb01b Mon Sep 17 00:00:00 2001 From: Eddie Schlafly Date: Thu, 14 Nov 2024 16:20:17 -0500 Subject: [PATCH 25/25] Get output filename correct. --- romancal/multiband_catalog/multiband_catalog_step.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/romancal/multiband_catalog/multiband_catalog_step.py b/romancal/multiband_catalog/multiband_catalog_step.py index 2d0079779..55c2b9fdb 100644 --- a/romancal/multiband_catalog/multiband_catalog_step.py +++ b/romancal/multiband_catalog/multiband_catalog_step.py @@ -68,6 +68,8 @@ def process(self, library): source_catalog_model.meta.filename = library.asn["products"][0]["name"] except (AttributeError, KeyError): source_catalog_model.meta.filename = "multiband_catalog" + if self.output_file is None: + self.output_file = source_catalog_model.meta.filename # TODO: sensible defaults # TODO: redefine in terms of intrinsic FWHM