From f22e4b9f68137bcca7f911c8c1496ff43cff0d6a Mon Sep 17 00:00:00 2001 From: mairan Date: Thu, 6 Jun 2024 20:56:05 -0400 Subject: [PATCH] RCAL-828: add PSF bits to source catalog. (#1243) Co-authored-by: zacharyburnett Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Brett Co-authored-by: D Davis <49163225+ddavis-stsci@users.noreply.github.com> --- CHANGES.rst | 7 +- docs/roman/source_catalog/arguments.rst | 21 +-- docs/roman/source_catalog/index.rst | 3 + docs/roman/source_catalog/main.rst | 41 ++++++ docs/roman/source_catalog/psf.rst | 29 ++++ docs/roman/source_detection/index.rst | 2 - romancal/lib/psf.py | 10 +- romancal/source_catalog/source_catalog.py | 136 +++++++++++++++++- .../source_catalog/source_catalog_step.py | 17 ++- .../tests/test_source_catalog.py | 131 ++++++++++++++--- 10 files changed, 353 insertions(+), 44 deletions(-) create mode 100644 docs/roman/source_catalog/psf.rst diff --git a/CHANGES.rst b/CHANGES.rst index 22f1c2c96..d92baa68b 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -14,11 +14,16 @@ general - Rename highlevelpipeline to mosaic pipeline [#1249] + +source_catalog +-------------- +- Add PSF photometry capability. [#1243] + dq_init ------- - - Refactor DQInitStep to use the RampModel method of creating ramps. [#1258] + 0.15.1 (2024-05-15) =================== diff --git a/docs/roman/source_catalog/arguments.rst b/docs/roman/source_catalog/arguments.rst index 26883c43e..5b11644aa 100644 --- a/docs/roman/source_catalog/arguments.rst +++ b/docs/roman/source_catalog/arguments.rst @@ -28,16 +28,19 @@ The ``source_catalog`` step uses the following user-settable arguments: energy value * ``--ci1_star_threshold``: A floating-point value of the threshold - compared to the concentration index calculated from aperture_ee1 - and aperture_ee2 that is used to determine whether a source is a - star. Sources must meet the criteria of both ci1_star_threshold and - ci2_star_threshold to be considered a star. + compared to the concentration index calculated from ``aperture_ee1`` + and ``aperture_ee2`` that is used to determine whether a source is a + star. Sources must meet the criteria of both ``ci1_star_threshold`` and + ``ci2_star_threshold`` to be considered a star. * ``--ci2_star_threshold``: A floating-point value of the threshold - compared to the concentration index calculated from aperture_ee2 - and aperture_ee3 that is used to determine whether a source is a - star. Sources must meet the criteria of both ci1_star_threshold and - ci2_star_threshold to be considered a star. + compared to the concentration index calculated from ``aperture_ee2`` + and ``aperture_ee3`` that is used to determine whether a source is a + star. Sources must meet the criteria of both ``ci1_star_threshold`` and + ``ci2_star_threshold`` to be considered a star. * ``--suffix``: A string value giving the file name suffix to use for - the output catalog file [default='cat'] + the output catalog file (default is ``'cat'``). + +* ``--fit_psf``: A boolean value indicating whether to perform + PSF photometry (default is ``True``) diff --git a/docs/roman/source_catalog/index.rst b/docs/roman/source_catalog/index.rst index 6c1221ee8..61a6cb085 100644 --- a/docs/roman/source_catalog/index.rst +++ b/docs/roman/source_catalog/index.rst @@ -9,5 +9,8 @@ Source Catalog main.rst arguments.rst + psf.rst .. automodapi:: romancal.source_catalog + +.. automodapi:: romancal.lib.psf diff --git a/docs/roman/source_catalog/main.rst b/docs/roman/source_catalog/main.rst index fea583d8e..7f0334c3e 100644 --- a/docs/roman/source_catalog/main.rst +++ b/docs/roman/source_catalog/main.rst @@ -60,6 +60,17 @@ Photometric errors are calculated from the resampled total-error array contained in the ``ERR`` (``model.err``) array. Note that this total-error array includes source Poisson noise. +PSF Fitting +----------- + +Star finding algorithms like `~photutils.detection.DAOStarFinder` provide +approximate stellar centroids. More precise centroids may be inferred by +fitting model PSFs to the observations. Setting the SourceDetectionStep's +option `fit_psf` to True will generate model Roman PSFs with +`WebbPSF `_, and fit +those models to each of the sources detected by +`~photutils.detection.DAOStarFinder`. More details are in :doc:`psf`. + Output Products --------------- @@ -170,6 +181,36 @@ columns (assuming the default encircled energies of 30, 50, and 70): | | minimal bounding box of the source | +------------------------+----------------------------------------------------+ + +If ``fit_psf=True``, the following columns will also be available: + ++------------------------+----------------------------------------------------+ +| Column | Description | ++========================+====================================================+ +| x_psf, x_psf_err | X pixel value of the source and its associated | +| | error as determined by PSF fitting | ++------------------------+----------------------------------------------------+ +| y_psf, y_psf_err | Y pixel value of the source and its associated | +| | error as determined by PSF fitting | ++------------------------+----------------------------------------------------+ +| flux_psf, flux_psf_err | Flux of the source and its associated error as | +| | determined by PSF fitting | ++------------------------+----------------------------------------------------+ +| flag_psf | DQ flag of the resulting PSF fitting. | +| | Possible values are [1]_: | +| | | +| | - 1 : one or more pixels in the fitting region | +| | were masked | +| | - 2 : the fit x and/or y position lies outside of | +| | the input data | +| | - 4 : the fit flux is less than or equal to zero | +| | - 8 : the fitter may not have converged | +| | - 16 : the fitter parameter covariance matrix was | +| | not returned | ++------------------------+----------------------------------------------------+ + +.. [1] See `PSFPhotometry `_ for more details. + Note that pixel coordinates are 0 indexed, matching the Python 0-based indexing. That means pixel coordinate ``0`` is the center of the first pixel. diff --git a/docs/roman/source_catalog/psf.rst b/docs/roman/source_catalog/psf.rst new file mode 100644 index 000000000..56f4af506 --- /dev/null +++ b/docs/roman/source_catalog/psf.rst @@ -0,0 +1,29 @@ +PSF Fitting +=========== + +A few PSF fitting utilities are included to interface between observations +within Roman datamodels and methods within dependencies that generate and +fit PSF models to observations. + +Create PSF models +----------------- + +`~romancal.lib.psf.create_gridded_psf_model` computes a gridded PSF model for +a given detector using `~webbpsf.gridded_library.CreatePSFLibrary` from +`WebbPSF `_. The defaults are chosen to +balance more accurate PSF models with the cost of increased runtime. For +further reading on the WebbPSF approach to ePSFs, see the WebbPSF docs on +`Using PSF Grids `_. + +Fit model PSFs to an ImageModel +------------------------------- + +Once PSF models are generated, you can fit those models to observations +in an ImageModel with `~romancal.lib.psf.fit_psf_to_image_model` using +`photutils `_. +By default the fits are done with `~photutils.psf.PSFPhotometry`, and +crowded fields may benefit from using `~photutils.psf.IterativePSFPhotometry`. +For neighboring sources that are near one another on the detector, grouping +the sources and fitting their PSFs simultaneously improves the fit quality. +Initial guesses for target centroids can be given or source +detection can be performed with, e.g., `~photutils.detection.DAOStarFinder`. diff --git a/docs/roman/source_detection/index.rst b/docs/roman/source_detection/index.rst index 31aa52bfd..2050333b9 100644 --- a/docs/roman/source_detection/index.rst +++ b/docs/roman/source_detection/index.rst @@ -12,5 +12,3 @@ Source Detection psf.rst .. automodapi:: romancal.source_detection - -.. automodapi:: romancal.lib.psf diff --git a/romancal/lib/psf.py b/romancal/lib/psf.py index f90ed921f..d95bcac6a 100644 --- a/romancal/lib/psf.py +++ b/romancal/lib/psf.py @@ -327,10 +327,12 @@ def fit_psf_to_image_model( # at some point in the future. if dq is None: - if image_model is not None: + if image_model is not None and isinstance(image_model, ImageModel): + # L2 images have a dq array mask = dq_to_boolean_mask(image_model, ignore_flags=ignore_flags) else: - mask = None + # L3 images + mask = image_model.weight == 0 else: mask = dq_to_boolean_mask(dq) @@ -352,9 +354,9 @@ def fit_psf_to_image_model( # outside the detector boundaries: init_centroid_in_range = ( (guesses["x_init"] > 0) - & (guesses["x_init"] < data.shape[0]) + & (guesses["x_init"] < data.shape[1]) & (guesses["y_init"] > 0) - & (guesses["y_init"] < data.shape[1]) + & (guesses["y_init"] < data.shape[0]) ) guesses = guesses[init_centroid_in_range] diff --git a/romancal/source_catalog/source_catalog.py b/romancal/source_catalog/source_catalog.py index 4437d360a..667283e08 100644 --- a/romancal/source_catalog/source_catalog.py +++ b/romancal/source_catalog/source_catalog.py @@ -4,6 +4,8 @@ import logging import warnings +from pathlib import Path +from typing import List import astropy.units as u import numpy as np @@ -22,6 +24,7 @@ from scipy.spatial import KDTree from romancal import __version__ as romancal_version +from romancal.lib import psf from ._wcs_helpers import pixel_scale_angle_at_skycoord @@ -84,6 +87,7 @@ def __init__( aperture_params, ci_star_thresholds, kernel_fwhm, + fit_psf, flux_unit="uJy", ): @@ -96,6 +100,7 @@ def __init__( self.aperture_params = aperture_params self.kernel_sigma = kernel_fwhm * gaussian_fwhm_to_sigma self.flux_unit = flux_unit + self.fit_psf = fit_psf self.sb_unit = "MJy/sr" self.l2_unit = "DN/s" @@ -113,6 +118,12 @@ def __init__( self.wcs = self.model.meta.wcs self.column_desc = {} self.meta = {} + if self.fit_psf: + # get the necessary columns for the PSF photometry table + # and its name mapping to the final catalog + self.psf_photometry_catalog_mapping = ( + self.get_psf_photometry_catalog_colnames_mapping() + ) @lazyproperty def _pixscale_angle(self): @@ -1007,10 +1018,12 @@ def colnames(self): """ The column name order for the final source catalog. """ - colnames = self.segment_colnames[0:4] + 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) return colnames def _update_metadata(self): @@ -1073,6 +1086,125 @@ def _split_skycoord(self, table): return table + @staticmethod + def get_psf_photometry_catalog_colnames_mapping() -> List: + """ + Set the mapping between the PSF photometry table column names + and the final catalog column names. + + Returns + ------- + List of dictionaries containing old column names, new column names, + and descriptions for PSF photometry catalog. + """ + + return [ + { + "old_name": "flags", + "new_name": "flag_psf", + "desc": "Data quality flags", + }, + { + "old_name": "x_fit", + "new_name": "x_psf", + "desc": "X coordinate as determined by PSF fitting", + }, + { + "old_name": "x_err", + "new_name": "x_psf_err", + "desc": "Error on X coordinate of PSF fitting", + }, + { + "old_name": "y_fit", + "new_name": "y_psf", + "desc": "Y coordinate as determined by PSF fitting", + }, + { + "old_name": "y_err", + "new_name": "y_psf_err", + "desc": "Error on Y coordinate of PSF fitting", + }, + { + "old_name": "flux_fit", + "new_name": "flux_psf", + "desc": "Source flux as determined by PSF photometry", + }, + { + "old_name": "flux_err", + "new_name": "flux_psf_err", + "desc": "Source flux error as determined by PSF photometry", + }, + ] + + @lazyproperty + def psf_photometry_colnames(self) -> List: + """ + Update and return column descriptions for PSF photometry results. + + This method updates the column descriptions with PSF photometry-related information + and returns a list of column names. + + Returns + ------- + List of column names for PSF photometry results. + """ + desc = {x["new_name"]: x["desc"] for x in self.psf_photometry_catalog_mapping} + + self.column_desc.update(desc) + + return list(desc.keys()) + + def do_psf_photometry(self) -> None: + """ + Perform PSF photometry by fitting PSF models to detected sources for refined astrometry. + + This method constructs a gridded PSF model based on instrument and detector information. + It then fits the PSF model to the image model's sources to improve astrometric precision. + + """ + log.info("Constructing a gridded PSF model.") + if hasattr(self.model.meta, "instrument"): + # ImageModel (L2 datamodel) + filt = self.model.meta.instrument.optical_element + detector = self.model.meta.instrument.detector.replace("WFI", "SCA") + else: + # MosaicModel (L3 datamodel) + filt = self.model.meta.basic.optical_element + detector = "SCA02" + # prefix of the temporary FITS file that will contain the gridded PSF model + path_prefix = "tmp" + gridded_psf_model, _ = psf.create_gridded_psf_model( + path_prefix=path_prefix, + filt=filt, + detector=detector, + overwrite=True, + logging_level="ERROR", + ) + + log.info("Fitting a PSF model to sources for improved astrometric precision.") + 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, + exclude_out_of_bounds=True, + ) + + # mapping between the columns of the PSF photometry table + # and the name that will be used in the final catalog + old_name_to_new_name_mapping = { + x["old_name"]: x["new_name"] for x in self.psf_photometry_catalog_mapping + } + + # append PSF results to the class instance with the proper column name + for old_name, new_name in old_name_to_new_name_mapping.items(): + setattr(self, new_name, psf_photometry_table[old_name]) + + # remove temporary file containing gridded_psf_model + filepath = Path().cwd().glob(f"{path_prefix}*{detector.lower()}*.fits") + for filename in filepath: + filename.unlink(missing_ok=True) + @lazyproperty def catalog(self): """ @@ -1086,6 +1218,8 @@ def catalog(self): self.set_segment_properties() self.set_aperture_properties() self.set_ci_properties() + if self.fit_psf: + self.do_psf_photometry() catalog = QTable() for column in self.colnames: diff --git a/romancal/source_catalog/source_catalog_step.py b/romancal/source_catalog/source_catalog_step.py index f1c811fb8..af72a5bff 100755 --- a/romancal/source_catalog/source_catalog_step.py +++ b/romancal/source_catalog/source_catalog_step.py @@ -42,6 +42,7 @@ class SourceCatalogStep(RomanStep): 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, input_model): @@ -58,15 +59,12 @@ def process(self, input_model): source_catalog_model = maker_utils.mk_datamodel(cat_model) for key in source_catalog_model.meta.keys(): - try: - if key == "optical_element": - value = model.meta.instrument[key] - else: - value = model.meta[key] - source_catalog_model.meta[key] = value - except KeyError: - pass - + value = ( + model.meta.instrument[key] + if key == "optical_element" + else model.meta[key] + ) + source_catalog_model.meta[key] = value aperture_ee = (self.aperture_ee1, self.aperture_ee2, self.aperture_ee3) refdata = ReferenceData(model, aperture_ee) aperture_params = refdata.aperture_params @@ -106,6 +104,7 @@ def process(self, input_model): aperture_params, ci_star_thresholds, self.kernel_fwhm, + self.fit_psf, ) # put the resulting catalog in the model diff --git a/romancal/source_catalog/tests/test_source_catalog.py b/romancal/source_catalog/tests/test_source_catalog.py index b95e7a061..cb58174e4 100644 --- a/romancal/source_catalog/tests/test_source_catalog.py +++ b/romancal/source_catalog/tests/test_source_catalog.py @@ -1,3 +1,5 @@ +import os + import astropy.units as u import numpy as np import pytest @@ -59,7 +61,7 @@ def make_test_image(): @pytest.fixture def mosaic_model(): - wfi_mosaic = mk_level3_mosaic() + wfi_mosaic = mk_level3_mosaic(shape=(101, 101)) model = MosaicModel(wfi_mosaic) data, err = make_test_image() units = u.MJy / u.sr @@ -67,12 +69,13 @@ def mosaic_model(): err <<= units model.data = data model.err = err + model.weight = 1.0 / err.value return model @pytest.fixture def image_model(): - wfi_image = mk_level2_image() + wfi_image = mk_level2_image(shape=(101, 101)) model = ImageModel(wfi_image) data, err = make_test_image() units = u.DN / u.s @@ -96,15 +99,19 @@ def image_model(): (50, 10, 0, False), ), ) -def test_l2_source_catalog(image_model, snr_threshold, npixels, nsources, save_results): - step = SourceCatalogStep( +def test_l2_source_catalog( + image_model, snr_threshold, npixels, nsources, save_results, tmp_path +): + os.chdir(tmp_path) + step = SourceCatalogStep() + result = step.call( + image_model, bkg_boxsize=50, kernel_fwhm=2.0, snr_threshold=snr_threshold, npixels=npixels, save_results=save_results, ) - result = step.run(image_model) cat = result.source_catalog assert isinstance(cat, Table) @@ -175,16 +182,19 @@ def test_l2_source_catalog(image_model, snr_threshold, npixels, nsources, save_r ), ) def test_l3_source_catalog( - mosaic_model, snr_threshold, npixels, nsources, save_results + mosaic_model, snr_threshold, npixels, nsources, save_results, tmp_path ): - step = SourceCatalogStep( + os.chdir(tmp_path) + step = SourceCatalogStep() + + result = step.call( + mosaic_model, bkg_boxsize=50, kernel_fwhm=2.0, snr_threshold=snr_threshold, npixels=npixels, save_results=save_results, ) - result = step.run(mosaic_model) cat = result.source_catalog assert isinstance(cat, Table) @@ -246,55 +256,61 @@ def test_background(mosaic_model): """ Test background fallback when Background2D fails. """ - step = SourceCatalogStep( + step = SourceCatalogStep() + result = step.call( + mosaic_model, bkg_boxsize=1000, kernel_fwhm=2.0, snr_threshold=3, npixels=25, + fit_psf=False, ) - result = step.run(mosaic_model) cat = result.source_catalog assert isinstance(cat, Table) -def test_l2_input_model_unchanged(image_model): +def test_l2_input_model_unchanged(image_model, tmp_path): """ Test that the input model data and error arrays are unchanged after processing by SourceCatalogStep. """ + os.chdir(tmp_path) original_data = image_model.data.copy() original_err = image_model.err.copy() - step = SourceCatalogStep( + step = SourceCatalogStep() + step.call( + image_model, snr_threshold=0.5, npixels=5, bkg_boxsize=50, kernel_fwhm=2.0, save_results=False, ) - step.run(image_model) assert_allclose(original_data, image_model.data, atol=5.0e-7) assert_allclose(original_err, image_model.err, atol=5.0e-7) -def test_l3_input_model_unchanged(mosaic_model): +def test_l3_input_model_unchanged(mosaic_model, tmp_path): """ Test that the input model data and error arrays are unchanged after processing by SourceCatalogStep. """ + os.chdir(tmp_path) original_data = mosaic_model.data.copy() original_err = mosaic_model.err.copy() - step = SourceCatalogStep( + step = SourceCatalogStep() + step.call( + mosaic_model, snr_threshold=0.5, npixels=5, bkg_boxsize=50, kernel_fwhm=2.0, save_results=False, ) - step.run(mosaic_model) assert_allclose(original_data, mosaic_model.data, atol=5.0e-7) assert_allclose(original_err, mosaic_model.err, atol=5.0e-7) @@ -323,7 +339,86 @@ def test_inputs(mosaic_model): aper_params = {} ci_thresh = 100.0 with pytest.raises(ValueError): - RomanSourceCatalog(np.ones((3, 3)), segm, cdata, aper_params, ci_thresh, 2.0) + RomanSourceCatalog( + np.ones((3, 3)), segm, cdata, aper_params, ci_thresh, 2.0, fit_psf=True + ) with pytest.raises(ValueError): - RomanSourceCatalog(mosaic_model, segm, cdata, aper_params, (1.0, 2.0, 3.0), 2.0) + RomanSourceCatalog( + mosaic_model, segm, cdata, aper_params, (1.0, 2.0, 3.0), 2.0, fit_psf=True + ) + + +def test_do_psf_photometry(tmp_path, image_model): + """ + Test that do_psf_photometry can recover mock sources and their position and photometry. + """ + os.chdir(tmp_path) + + # get column names mapping for PSF photometry + psf_colnames_mapping = ( + RomanSourceCatalog.get_psf_photometry_catalog_colnames_mapping() + ) + psf_colnames = [ + x.get("new_name") + for x in psf_colnames_mapping + if x.get("old_name") in ["x_fit", "y_fit", "flux_fit"] + ] + + step = SourceCatalogStep() + result = step.call( + image_model, + bkg_boxsize=20, + kernel_fwhm=2.0, + snr_threshold=3, + npixels=10, + save_results=False, + ) + cat = result.source_catalog + + assert isinstance(cat, Table) + + # check the number of sources that have been detected + assert len(cat) == 7 + # check that all sources have both position and flux determined (ignore errors/flags) + assert all(len(cat[x]) and cat[x] is not [None, np.nan] for x in psf_colnames) + + +@pytest.mark.parametrize("fit_psf", [True, False]) +def test_do_psf_photometry_column_names(tmp_path, image_model, fit_psf): + """ + Test that fit_psf will determine whether the PSF + photometry columns are added to the final catalog or not. + """ + os.chdir(tmp_path) + + # get column names mapping for PSF photometry + psf_colnames_mapping = ( + RomanSourceCatalog.get_psf_photometry_catalog_colnames_mapping() + ) + + step = SourceCatalogStep() + result = step.call( + image_model, + bkg_boxsize=20, + kernel_fwhm=2.0, + snr_threshold=3, + npixels=10, + save_results=False, + fit_psf=fit_psf, + ) + cat = result.source_catalog + + assert isinstance(cat, Table) + + # check if the PSF photometry column names are present or not based on fit_psf value + psf_colnames_present = all( + x.get("new_name") in cat.colnames for x in psf_colnames_mapping + ) + psf_colnames_not_present = all( + x.get("new_name") not in cat.colnames for x in psf_colnames_mapping + ) + + assert (fit_psf and psf_colnames_present) or ( + not fit_psf and psf_colnames_not_present + )