From 2c5431388f32877d1d8031bfbec6709bbacac12d Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Thu, 24 Jul 2025 21:11:22 -0700 Subject: [PATCH 01/23] Check if column is masked before consolidating --- .../_consolidate_injected_coadd_catalogs.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/python/lsst/source/injection/utils/_consolidate_injected_coadd_catalogs.py b/python/lsst/source/injection/utils/_consolidate_injected_coadd_catalogs.py index bb9f00c..2ed7fd2 100644 --- a/python/lsst/source/injection/utils/_consolidate_injected_coadd_catalogs.py +++ b/python/lsst/source/injection/utils/_consolidate_injected_coadd_catalogs.py @@ -33,7 +33,7 @@ import astropy.units as u import numpy as np from astropy.table import Table, join, vstack -from astropy.table.column import MaskedColumn +from astropy.table.column import Column, MaskedColumn from smatch.matcher import Matcher # type: ignore [import-not-found] from lsst.daf.butler import DatasetProvenance @@ -546,15 +546,22 @@ def make_multiband_catalog( # One component is sufficient to say it was injected injection_flag[idx_old] &= not injected + # Convert any masked columns to regular nan-valued arrays if is_multicomp: column = values_comp[f"{band}_{prefix_comp}flux"] - column.mask[idx_old] = False - column[idx_old] = flux_comp + if isinstance(column, MaskedColumn): + column.mask[idx_old] = False + column[idx_old] = Column(column) + else: + column[idx_old] = flux_comp for column_comp in columns_comp: - column = f"{band}_{prefix_comp}{column_comp}" - values_comp[column].mask[idx_old] = False - values_comp[column][idx_old] = row[column_comp] + column = values_comp[f"{band}_{prefix_comp}{column_comp}"] + if isinstance(column, MaskedColumn): + column.mask[idx_old] = False + column[idx_old] = Column(column) + else: + column[idx_old] = row[column_comp] if is_multicomp: mag = (flux * u.nJy).to(u.ABmag).value From 0db5e58050545492bcced1cc8e171ac855999ed4 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Thu, 24 Jul 2025 21:11:48 -0700 Subject: [PATCH 02/23] Add parameters to post-injection pipeline --- pipelines/match_injected_tract_catalog.yaml | 72 ++++++++++++++------- 1 file changed, 47 insertions(+), 25 deletions(-) diff --git a/pipelines/match_injected_tract_catalog.yaml b/pipelines/match_injected_tract_catalog.yaml index 50953fb..3382f88 100644 --- a/pipelines/match_injected_tract_catalog.yaml +++ b/pipelines/match_injected_tract_catalog.yaml @@ -1,4 +1,8 @@ description: Match tract-level injected coadd catalogs to injected object table datasets. +parameters: + bands_injected: ["u", "g", "r", "i", "z", "y"] + flux_algorithms: ["cModel", "gaap1p0", "gaap3p0", "kron", "psf", "sersic"] + source_type_column: "source_type" tasks: consolidateInjectedCatalogs: class: lsst.source.injection.utils.ConsolidateInjectedCatalogsTask @@ -29,30 +33,48 @@ tasks: connections.name_input_cat_ref: injected_deep_coadd_predetection_catalog_tract connections.name_input_cat_target: injected_objectTable_tract column_matched_prefix_ref: "ref_" - # TODO: Remove as part of DM-44139 - columns_ref_mag_to_nJy: { - "u_mag": "u_flux", "g_mag": "g_flux", "r_mag": "r_flux", - "i_mag": "i_flux", "z_mag": "z_flux", "y_mag": "y_flux", - } - columns_ref_copy: [ - "injected_id", "source_type", - "u_mag", "g_mag", "r_mag", "i_mag", "z_mag", "y_mag", - "u_injection_flag", "g_injection_flag", "r_injection_flag", - "i_injection_flag", "z_injection_flag", "y_injection_flag", - "injection_flag", "injected_isPatchInner", "injected_isTractInner", "injected_isPrimary", - ] columns_target_coord_err: ["coord_raErr", "coord_decErr"] - columns_target_copy: [ - "objectId", - "u_psfFlux", "g_psfFlux", "r_psfFlux", "i_psfFlux", "z_psfFlux", "y_psfFlux", - "u_psfFluxErr", "g_psfFluxErr", "r_psfFluxErr", "i_psfFluxErr", "z_psfFluxErr", "y_psfFluxErr", - "u_psfFlux_flag", "g_psfFlux_flag", "r_psfFlux_flag", - "i_psfFlux_flag", "z_psfFlux_flag", "y_psfFlux_flag", - "u_sersicFlux", "g_sersicFlux", "r_sersicFlux", "i_sersicFlux", "z_sersicFlux", "y_sersicFlux", - "u_sersicFluxErr", "g_sersicFluxErr", "r_sersicFluxErr", - "i_sersicFluxErr", "z_sersicFluxErr", "y_sersicFluxErr", - "sersic_ra", "sersic_dec", - "patch", "detect_isDeblendedSource", "detect_isPatchInner", "detect_isPrimary", "merge_peak_sky", - "refExtendedness", "refSizeExtendedness", - ] include_unmatched: true + python: | + columns_ref_mag_to_nJy = {} + columns_ref_copy = [ + "injected_id", parameters.source_type_column, "injection_flag", + "injected_isPatchInner", "injected_isTractInner", "injected_isPrimary", + ] + columns_target_copy = [ + "objectId", "patch", "detect_isDeblendedSource", "detect_isPatchInner", "detect_isPrimary", + "merge_peak_sky", "refExtendedness", "refSizeExtendedness", + ] + bands = parameters.bands_injected + for band in bands: + columns_ref_mag_to_nJy[f"{band}_mag"] = f"{band}_flux" + for band in bands: + for suffix in ("mag", "injection_flag"): + columns_ref_copy.append(f"{band}_{suffix}") + done_gaap = False + for algo in parameters.flux_algorithms: + is_gaap = algo.startswith("gaap") + if not done_gaap and is_gaap: + columns_target_copy.append(f"{band}_gaapFlux_flag") + done_gaap = True + suffixes = ["Flux", "FluxErr"] + if is_gaap: + suffixes.append("Flux_flag_bigPsf") + elif algo == "cModel": + suffixes.append("_flag") + elif algo != "sersic": + suffixes.append("Flux_flag") + for suffix in suffixes: + columns_target_copy.append(f"{band}_{algo}{suffix}") + if "sersic" in parameters.flux_algorithms: + columns_target_copy.extend([ + "sersic_x", "sersic_y", "sersic_xErr", "sersic_yErr", + "sersic_ra", "sersic_dec", "sersic_raErr", "sersic_decErr", + "sersic_index", "sersic_indexErr", "sersic_rho", "sersic_rhoErr", + "sersic_reff_x", "sersic_reff_xErr", "sersic_reff_y", "sersic_reff_yErr", + "sersic_unknown_flag", "sersic_no_data_flag", + ]) + # TODO: Remove as part of DM-44139 + config.columns_ref_mag_to_nJy = columns_ref_mag_to_nJy + config.columns_ref_copy = columns_ref_copy + config.columns_target_copy = columns_target_copy From b4d71b46ce3c02fac7502722cd2c87c5868b02bd Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Fri, 10 Oct 2025 10:56:51 -0700 Subject: [PATCH 03/23] Simplify matching with DM-51163 --- pipelines/match_injected_tract_catalog.yaml | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/pipelines/match_injected_tract_catalog.yaml b/pipelines/match_injected_tract_catalog.yaml index 3382f88..727f72a 100644 --- a/pipelines/match_injected_tract_catalog.yaml +++ b/pipelines/match_injected_tract_catalog.yaml @@ -23,19 +23,18 @@ tasks: match_tract_catalog.columns_target_select_false: ["merge_peak_sky"] match_tract_catalog.match_n_finite_min: 2 match_tract_catalog.order_ascending: true + output_matched_catalog: true + + diff_matched_catalog.column_matched_prefix_ref: "ref_" + diff_matched_catalog.columns_target_coord_err: ["coord_raErr", "coord_decErr"] + diff_matched_catalog.columns_ref_copy: columns_ref_copy + diff_matched_catalog.columns_target_copy: columns_target_copy + diff_matched_catalog.include_unmatched: true + python: | from lsst.pipe.tasks.match_tract_catalog_probabilistic import MatchTractCatalogProbabilisticTask config.match_tract_catalog.retarget(MatchTractCatalogProbabilisticTask) - compareObjectToInjected: - class: lsst.pipe.tasks.diff_matched_tract_catalog.DiffMatchedTractCatalogTask - config: - connections.name_input_cat_ref: injected_deep_coadd_predetection_catalog_tract - connections.name_input_cat_target: injected_objectTable_tract - column_matched_prefix_ref: "ref_" - columns_target_coord_err: ["coord_raErr", "coord_decErr"] - include_unmatched: true - python: | columns_ref_mag_to_nJy = {} columns_ref_copy = [ "injected_id", parameters.source_type_column, "injection_flag", @@ -75,6 +74,4 @@ tasks: "sersic_unknown_flag", "sersic_no_data_flag", ]) # TODO: Remove as part of DM-44139 - config.columns_ref_mag_to_nJy = columns_ref_mag_to_nJy - config.columns_ref_copy = columns_ref_copy - config.columns_target_copy = columns_target_copy + config.diff_matched_catalog.columns_ref_mag_to_nJy = columns_ref_mag_to_nJy From 77a211d6dc0c1b3bc42f01806ab60f696e485d4f Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Wed, 9 Sep 2026 11:18:03 -0700 Subject: [PATCH 04/23] Update black/python version in pre-commit --- .pre-commit-config.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 322f32e..c3b40a8 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,14 +6,14 @@ repos: - id: end-of-file-fixer - id: trailing-whitespace - repo: https://github.com/psf/black - rev: 24.4.2 + rev: 25.12.0 hooks: - id: black # It is recommended to specify the latest version of Python # supported by your project here, or alternatively use # pre-commit's default_language_version, see # https://pre-commit.com/#top_level-default_language_version - language_version: python3.11 + language_version: python3.13 - repo: https://github.com/PyCQA/isort rev: 5.13.2 hooks: From d3ff258dd3739884e781d8237d0f22803fd414cb Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Tue, 7 Jul 2026 20:03:12 -0700 Subject: [PATCH 05/23] Move PSF image computation to _get_galsim_psf This is mainly to improve the readibility of a very long function that will need more specialization for cell coadds. --- python/lsst/source/injection/inject_engine.py | 68 ++++++++++++------- 1 file changed, 42 insertions(+), 26 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index 63e8a45..a465228 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -597,6 +597,34 @@ def add_noise_to_galsim_image( return image_template, var_template +def _get_galsim_psf(*, psf, pixel_coords, bbox, calib_flux_radius, galsim_wcs, sky_coords=None): + """Compute the galsim PSF at given coordinates within a bbox.""" + psf_array = None + error_args = [] + try: + psf_array = psf.computeKernelImage(pixel_coords).array + except InvalidParameterError: + # Try mapping to nearest point contained in bbox. + contained_point = Point2D( + np.clip(pixel_coords.x, bbox.minX, bbox.maxX), np.clip(pixel_coords.y, bbox.minY, bbox.maxY) + ) + if pixel_coords == contained_point: # no difference, so skip immediately + error_args = ["Cannot compute PSF for object at %s; flagging and skipping.", sky_coords] + # Otherwise, try again with new point. + try: + psf_array = psf.computeKernelImage(contained_point).array + except InvalidParameterError: + error_args = ["Cannot compute PSF for object at %s; flagging and skipping.", sky_coords] + + if psf_array is not None: + # Compute the aperture corrected PSF interpolated image. + aperture_correction = psf.computeApertureFlux(calib_flux_radius, psf.getAveragePosition()) + psf_array /= aperture_correction + psf_array = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) + + return psf_array, error_args + + def inject_galsim_objects_into_exposure( exposure: ExposureF, objects: Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None], @@ -702,32 +730,20 @@ def inject_galsim_objects_into_exposure( if galsim_pixel_scale < pixel_scale / 2 or galsim_pixel_scale > pixel_scale * 2: continue - # Compute the PSF at the object location. - try: - psf_array = psf.computeKernelImage(pixel_coords).array - except InvalidParameterError: - # Try mapping to nearest point contained in bbox. - contained_point = Point2D( - np.clip(pixel_coords.x, bbox.minX, bbox.maxX), np.clip(pixel_coords.y, bbox.minY, bbox.maxY) - ) - if pixel_coords == contained_point: # no difference, so skip immediately - psf_compute_errors[i] = True - if logger: - logger.debug("Cannot compute PSF for object at %s; flagging and skipping.", sky_coords) - continue - # Otherwise, try again with new point. - try: - psf_array = psf.computeKernelImage(contained_point).array - except InvalidParameterError: - psf_compute_errors[i] = True - if logger: - logger.debug("Cannot compute PSF for object at %s; flagging and skipping.", sky_coords) - continue - - # Compute the aperture corrected PSF interpolated image. - aperture_correction = psf.computeApertureFlux(calib_flux_radius, psf.getAveragePosition()) - psf_array /= aperture_correction - galsim_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) + # Get the PSF at the centroid of the object + galsim_psf_centroid, aperture_correction, psf_error_args = _get_galsim_psf( + psf=psf, + pixel_coords=pixel_coords, + bbox=bbox, + calib_flux_radius=calib_flux_radius, + galsim_wcs=galsim_wcs, + sky_coords=sky_coords, + ) + if galsim_psf is None: + psf_compute_errors[i] = True + if logger: + logger.debug(*psf_error_args) + continue # Convolve the object with the PSF and generate draw size. conv = galsim.Convolve(object, galsim_psf) From a12509dcc32637cabdbaba47cb581372ff6d97bf Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Wed, 26 Aug 2026 21:05:18 -0700 Subject: [PATCH 06/23] Move object injection into separate method This shortens the very long inject_galsim_objects_into_exposure method in advance of it growing again to support cell coadds, albeit at the expense of adding a private method with many (keyword) arguments. --- python/lsst/source/injection/inject_engine.py | 151 +++++++++++------- 1 file changed, 97 insertions(+), 54 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index a465228..34cba63 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -625,6 +625,84 @@ def _get_galsim_psf(*, psf, pixel_coords, bbox, calib_flux_radius, galsim_wcs, s return psf_array, error_args +def _inject_galsim_object_into_bounds( + *, + convolved_object, + position_d, + position_i, + object_common_bounds, + full_bounds, + galsim_image: galsim.Image, + galsim_variance: galsim.Image, + galsim_wcs, + gain_map, + exposure, + inject_variance: bool, + add_noise: bool, + noise_seed: int, + is_coadd: bool, + injection_core_size, + mask_plane_name, + mask_plane_core_name, + logger, +): + common_image = galsim_image[object_common_bounds] + common_variance = galsim_variance[object_common_bounds] + + offset = position_d - object_common_bounds.true_center + # Note, for preliminary_visit_image injection, pixel is already + # part of the PSF and for coadd injection, it's incorrect to + # include the output pixel. + # So for both cases, we draw using method='no_pixel'. + image_template = common_image.copy() + image_template = convolved_object.drawImage( + image_template, add_to_image=False, offset=offset, wcs=galsim_wcs, method="no_pixel" + ) + + # The amount of additional variance to inject, + # if inject_variance is true. + var_template = image_template.copy() + var_template *= gain_map[object_common_bounds] + + if add_noise: + image_template, var_template = add_noise_to_galsim_image( + image_template, + var_template, + is_coadd, + gain_map, + object_common_bounds, + noise_seed, + logger, + ) + + common_image += image_template + if inject_variance: + common_variance += var_template + + # Increment the seed so different noise is generated for different + # objects. + noise_seed += 1 + + common_box = Box2I( + Point2I(object_common_bounds.xmin, object_common_bounds.ymin), + Point2I(object_common_bounds.xmax, object_common_bounds.ymax), + ) + bitvalue = exposure.mask.getPlaneBitMask(mask_plane_name) + exposure[common_box].mask.array |= bitvalue + # Add a configurable pixel mask centered on the object. The mask must be + # large enough to always identify the core/peak of the injected + # source, but small enough that it rarely overlaps real sources. + sub_bounds_core = galsim.BoundsI(position_i).withBorder(injection_core_size // 2) + object_common_bounds_core = full_bounds & sub_bounds_core + if object_common_bounds_core.area() > 0: + common_box_core = Box2I( + Point2I(object_common_bounds_core.xmin, object_common_bounds_core.ymin), + Point2I(object_common_bounds_core.xmax, object_common_bounds_core.ymax), + ) + bitvalue_core = exposure.mask.getPlaneBitMask(mask_plane_core_name) + exposure[common_box_core].mask.array |= bitvalue_core + + def inject_galsim_objects_into_exposure( exposure: ExposureF, objects: Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None], @@ -776,18 +854,26 @@ def inject_galsim_objects_into_exposure( # Inject the source if there is any overlap. if object_common_bounds.area() > 0: - common_image = galsim_image[object_common_bounds] - common_variance = galsim_variance[object_common_bounds] - - offset = posd - object_common_bounds.true_center - # Note, for preliminary_visit_image injection, pixel is already - # part of the PSF and for coadd injection, it's incorrect to - # include the output pixel. - # So for both cases, we draw using method='no_pixel'. - image_template = common_image.copy() try: - image_template = conv.drawImage( - image_template, add_to_image=False, offset=offset, wcs=galsim_wcs, method="no_pixel" + _inject_galsim_object_into_bounds( + convolved_object=conv, + position_d=posd, + position_i=posi, + object_common_bounds=object_common_bounds, + full_bounds=full_bounds, + galsim_image=galsim_image, + galsim_variance=galsim_variance, + galsim_wcs=galsim_wcs, + gain_map=gain_map, + exposure=exposure, + inject_variance=inject_variance, + add_noise=add_noise, + noise_seed=noise_seed, + is_coadd=is_coadd, + injection_core_size=injection_core_size, + mask_plane_name=mask_plane_name, + mask_plane_core_name=mask_plane_core_name, + logger=logger, ) except GalSimFFTSizeError as err: fft_size_errors[i] = True @@ -798,49 +884,6 @@ def inject_galsim_objects_into_exposure( err, ) continue - - # The amount of additional variance to inject, - # if inject_variance is true. - var_template = image_template.copy() - var_template *= gain_map[object_common_bounds] - - if add_noise: - image_template, var_template = add_noise_to_galsim_image( - image_template, - var_template, - is_coadd, - gain_map, - object_common_bounds, - noise_seed, - logger, - ) - - common_image += image_template - if inject_variance: - common_variance += var_template - - # Increment the seed so different noise is generated for different - # objects. - noise_seed += 1 - - common_box = Box2I( - Point2I(object_common_bounds.xmin, object_common_bounds.ymin), - Point2I(object_common_bounds.xmax, object_common_bounds.ymax), - ) - bitvalue = exposure.mask.getPlaneBitMask(mask_plane_name) - exposure[common_box].mask.array |= bitvalue - # Add a 3 x 3 pixel mask centered on the object. The mask must be - # large enough to always identify the core/peak of the injected - # source, but small enough that it rarely overlaps real sources. - sub_bounds_core = galsim.BoundsI(posi).withBorder(injection_core_size // 2) - object_common_bounds_core = full_bounds & sub_bounds_core - if object_common_bounds_core.area() > 0: - common_box_core = Box2I( - Point2I(object_common_bounds_core.xmin, object_common_bounds_core.ymin), - Point2I(object_common_bounds_core.xmax, object_common_bounds_core.ymax), - ) - bitvalue_core = exposure.mask.getPlaneBitMask(mask_plane_core_name) - exposure[common_box_core].mask.array |= bitvalue_core else: if logger: logger.debug("No area overlap for object at %s; flagging and skipping.", sky_coords) From f99e9ed6edcfbf37830c0e1ec03b523ea8cab4e6 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Wed, 26 Aug 2026 21:16:47 -0700 Subject: [PATCH 07/23] Implement injecting into cell coadds --- python/lsst/source/injection/inject_engine.py | 177 ++++++++++++------ 1 file changed, 124 insertions(+), 53 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index 34cba63..f281411 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -29,6 +29,7 @@ from typing import Any import galsim +import galsim.errors as galsim_errors import numpy as np import numpy.ma as ma from astropy.io import fits @@ -38,8 +39,13 @@ from lsst.afw.geom import SkyWcs from lsst.afw.image import ExposureF, PhotoCalib from lsst.geom import Box2I, Point2D, Point2I, SpherePoint, arcseconds, degrees +from lsst.images import BoundsError +from lsst.images.cells import CellCoadd, CellIJ from lsst.pex.exceptions import InvalidParameterError, LogicError +# This is off by default but we want it on since we catch the error +galsim_errors.raise_fft_size_error = True + def get_object_data(source_data: dict[str, Any], object_class: galsim.GSObject) -> dict[str, Any]: """Assemble a dictionary of allowed keyword arguments and their @@ -618,6 +624,7 @@ def _get_galsim_psf(*, psf, pixel_coords, bbox, calib_flux_radius, galsim_wcs, s if psf_array is not None: # Compute the aperture corrected PSF interpolated image. + # TODO: Review whether we want to at least make this optional aperture_correction = psf.computeApertureFlux(calib_flux_radius, psf.getAveragePosition()) psf_array /= aperture_correction psf_array = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) @@ -679,10 +686,6 @@ def _inject_galsim_object_into_bounds( if inject_variance: common_variance += var_template - # Increment the seed so different noise is generated for different - # objects. - noise_seed += 1 - common_box = Box2I( Point2I(object_common_bounds.xmin, object_common_bounds.ymin), Point2I(object_common_bounds.xmax, object_common_bounds.ymax), @@ -704,7 +707,7 @@ def _inject_galsim_object_into_bounds( def inject_galsim_objects_into_exposure( - exposure: ExposureF, + exposure: ExposureF | CellCoadd, objects: Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None], mask_plane_name: str = "INJECTED", calib_flux_radius: float = 12.0, @@ -719,7 +722,7 @@ def inject_galsim_objects_into_exposure( Parameters ---------- - exposure : `lsst.afw.image.ExposureF` + exposure : `lsst.afw.image.ExposureF` or `lsst.images.cells.CellCoadd` The exposure to inject synthetic sources into. objects : `Generator` [`tuple`, None, None] An iterator of tuples that contains (or generates) locations and object @@ -763,25 +766,62 @@ def inject_galsim_objects_into_exposure( psf_compute_errors : `list` [`bool`] Boolean flags indicating whether a PSF computation error was raised. """ + is_cell = isinstance(exposure, CellCoadd) + if not is_cell and not isinstance(exposure, ExposureF): + raise TypeError(f"Unsupported {type(exposure)=}") + if is_cell: + cell_coadd = exposure + psf = exposure.psf + exposure = cell_coadd.to_legacy() + is_coadd = True + else: + is_coadd = exposure.info.getCoaddInputs() is not None + exposure.mask.addMaskPlane(mask_plane_name) mask_plane_core_name = mask_plane_name + "_CORE" exposure.mask.addMaskPlane(mask_plane_core_name) + wcs = exposure.getWcs() + bbox = exposure.getBBox() if logger: logger.info( "Adding %s and %s mask planes to the exposure.", mask_plane_name, mask_plane_core_name, ) - psf = exposure.getPsf() - wcs = exposure.getWcs() - bbox = exposure.getBBox() full_bounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY) galsim_image = galsim.Image(exposure.image.array, bounds=full_bounds) galsim_variance = galsim.Image(exposure.variance.array, bounds=full_bounds) pixel_scale = wcs.getPixelScale(bbox.getCenter()).asArcseconds() + if is_cell: + all_bounds = [] + cell_grid = cell_coadd.grid + n_i, n_j = cell_grid.grid_size.as_tuple() + + fallback_psf = None + + # Iterate over all cells and find the ones with good PSFs + for cell_i in range(n_i): + for cell_j in range(n_j): + cell_ij = CellIJ(cell_i, cell_j) + bbox_cell = cell_grid.bbox_of(cell_ij) + cen_bbox = bbox_cell.to_legacy().getCenter() + cen_x, cen_y = cen_bbox + try: + psf_array = cell_coadd.psf.compute_kernel_image(x=cen_x, y=cen_y) + all_bounds.append((bbox_cell, (cell_ij, cen_x, cen_y))) + if fallback_psf is None: + mat = wcs.linearizePixelToSky(cen_bbox).getMatrix() + galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) + fallback_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) + except BoundsError: + # Assume this cell has a bad PSF/data and give up + pass + else: + all_bounds = [(full_bounds, None)] + + # TODO: Refactor to not take exposure? gain_map = get_gain_map(exposure, bad_mask_names, logger) - is_coadd = exposure.info.getCoaddInputs() is not None draw_sizes: list[int] = [] common_bounds: list[galsim.BoundsI] = [] @@ -809,7 +849,7 @@ def inject_galsim_objects_into_exposure( continue # Get the PSF at the centroid of the object - galsim_psf_centroid, aperture_correction, psf_error_args = _get_galsim_psf( + galsim_psf_centroid, psf_error_args = _get_galsim_psf( psf=psf, pixel_coords=pixel_coords, bbox=bbox, @@ -817,14 +857,27 @@ def inject_galsim_objects_into_exposure( galsim_wcs=galsim_wcs, sky_coords=sky_coords, ) - if galsim_psf is None: - psf_compute_errors[i] = True - if logger: - logger.debug(*psf_error_args) - continue + # Get the fallback PSF (if available) to compute the injection box. + # Even if the PSF at the centroid of an object is not valid, nearby + # cells might be fine. For example, cells with saturated stars can + # end up with no visits below the masked pixel threshold and no PSF. + # + # Note that _get_galsim_psf already implements a fallback of getting + # the PSF at the nearest edge of the bbox for objects whose centroids + # land outside the exposure's bbox. Non-cell coadds and visits could + # also do something similar by searching for the nearest valid PSF, + # but that's not as trivial as picking an "average" cell. + if galsim_psf_centroid is None: + if is_cell: + galsim_psf_centroid, psf_error_args = fallback_psf, [] + else: + psf_compute_errors[i] = True + if logger: + logger.debug(*psf_error_args) + continue # Convolve the object with the PSF and generate draw size. - conv = galsim.Convolve(object, galsim_psf) + conv = galsim.Convolve(object, galsim_psf_centroid) if draw_size == 0: draw_size = conv.getGoodImageSize(galsim_wcs.minLinearScale()) # type: ignore injection_draw_size = int(draw_size) @@ -849,43 +902,61 @@ def inject_galsim_objects_into_exposure( ) injection_core_size = injection_draw_size sub_bounds = galsim.BoundsI(posi).withBorder(injection_draw_size // 2) - object_common_bounds = full_bounds & sub_bounds - common_bounds[i] = object_common_bounds # type: ignore - - # Inject the source if there is any overlap. - if object_common_bounds.area() > 0: - try: - _inject_galsim_object_into_bounds( - convolved_object=conv, - position_d=posd, - position_i=posi, - object_common_bounds=object_common_bounds, - full_bounds=full_bounds, - galsim_image=galsim_image, - galsim_variance=galsim_variance, - galsim_wcs=galsim_wcs, - gain_map=gain_map, - exposure=exposure, - inject_variance=inject_variance, - add_noise=add_noise, - noise_seed=noise_seed, - is_coadd=is_coadd, - injection_core_size=injection_core_size, - mask_plane_name=mask_plane_name, - mask_plane_core_name=mask_plane_core_name, - logger=logger, - ) - except GalSimFFTSizeError as err: - fft_size_errors[i] = True - if logger: - logger.debug( - "GalSimFFTSizeError raised for object at %s; flagging and skipping.\n%s", - sky_coords, - err, + + # These bounds may not be the same as what would be derived from + # per-cell PSFs for large objects, but re-calculating them for every + # cell is not worthwhile + if is_cell: + object_common_bounds = full_bounds & sub_bounds + if object_common_bounds.area() > 0: + common_bounds[i] = object_common_bounds # type: ignore + + area_to_inject = 0 + area_injected = 0 + bounds_fft_size_errors = {} + for idx_bound, (injection_bounds, bounds_info) in enumerate(all_bounds): + object_common_bounds = injection_bounds & sub_bounds + + # Inject the source if there is any overlap. + if (common_area := object_common_bounds.area()) > 0: + area_to_inject += common_area + if not is_cell: + common_bounds[i] = object_common_bounds # type: ignore + try: + _inject_galsim_object_into_bounds( + convolved_object=conv, + position_d=posd, + position_i=posi, + object_common_bounds=object_common_bounds, + full_bounds=full_bounds, + galsim_image=galsim_image, + galsim_variance=galsim_variance, + galsim_wcs=galsim_wcs, + gain_map=gain_map, + exposure=exposure, + inject_variance=inject_variance, + add_noise=add_noise, + noise_seed=noise_seed, + is_coadd=is_coadd, + injection_core_size=injection_core_size, + mask_plane_name=mask_plane_name, + mask_plane_core_name=mask_plane_core_name, + logger=logger, ) - continue - else: + area_injected += 0 + except GalSimFFTSizeError as err: + bounds_fft_size_errors[idx_bound] = err.size + if logger and (area_to_inject == 0): + logger.debug("No area overlap for object at %s; flagging and skipping.", sky_coords) + if bounds_fft_size_errors: + fft_size_errors[i] = True if logger: - logger.debug("No area overlap for object at %s; flagging and skipping.", sky_coords) + logger.debug( + "GalSimFFTSizeError raised for object at index=%d and coords %s;" + " flagging and skipping.\nbounds_fft_size_errors=%s", + i, + sky_coords, + bounds_fft_size_errors, + ) return draw_sizes, common_bounds, fft_size_errors, psf_compute_errors From 8276ffa4936fc9a4b66e920d32f649cf385f784e Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Thu, 27 Aug 2026 17:22:10 -0700 Subject: [PATCH 08/23] Add make_test_cell_coadd to convert legacy exposures --- .../lsst/source/injection/utils/test_utils.py | 96 ++++++++++++++++++- 1 file changed, 94 insertions(+), 2 deletions(-) diff --git a/python/lsst/source/injection/utils/test_utils.py b/python/lsst/source/injection/utils/test_utils.py index d398a3e..97c0df1 100644 --- a/python/lsst/source/injection/utils/test_utils.py +++ b/python/lsst/source/injection/utils/test_utils.py @@ -19,12 +19,22 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . - +import astropy.units as u import numpy as np from lsst.afw.geom import makeCdMatrix, makeSkyWcs -from lsst.afw.image import makePhotoCalibFromCalibZeroPoint +from lsst.afw.image import ExposureF, makePhotoCalibFromCalibZeroPoint from lsst.geom import Box2I, Extent2I, Point2D, Point2I, SpherePoint, degrees +from lsst.images import YX, Box, Image, Mask, SkyProjection, TractFrame +from lsst.images.cells import ( + CellCoadd, + CellGrid, + CellGridBounds, + CellIJ, + CellPointSpreadFunction, + CoaddProvenance, + PatchDefinition, +) from lsst.ip.isr.isrTask import IsrTask from lsst.meas.algorithms.testUtils import plantSources from lsst.pipe.base import Pipeline @@ -34,6 +44,88 @@ from lsst.source.injection import generate_injection_catalog +def make_test_cell_coadd( + exposure: ExposureF, + cell_shape: YX, + psf_shape: tuple[int, int], + missing: set[CellIJ] | None = None, + image_unit: u.Unit = u.nJy, + **kwargs, +) -> CellCoadd: + """Make a cell coadd out of an exposure. + + Parameters + ---------- + exposure : `lsst.afw.image.ExposureF` + An exposure with a simple PSF. + cell_shape : `lsst.images.YX` + The shape of each cell. + psf_shape : `tuple[int, int]` + The shape of the PSF image array. + missing : `set[CellIJ]` or `None` + The set of cells that have invalid PSFs, which will be set to nan. + image_unit : `astropy.units.Unit` + The unit of the exposure's image plane. + kwargs + Additional keyword arguments to pass to the CellCoadd constructor. + + Returns + ------- + cell_coadd : `lsst.images.cells.CellCoadd` + The exposure as a CellCoadd. + """ + if missing is None: + missing = set() + + exposure_psf = exposure.psf + box = Box.from_legacy(exposure.getBBox()) + cell_grid = CellGrid(bbox=box, cell_shape=cell_shape) + cell_grid_bounds = CellGridBounds(grid=cell_grid, bbox=box, missing=missing) + + # Make the cell PSF grid by evaluating the PSF at the center of each cell + cell_psfs = [] + for cell_i in range(cell_grid.grid_size.i): + row_psfs = [] + for cell_j in range(cell_grid.grid_size.j): + if (cell_ij := CellIJ(cell_i, cell_j)) in missing: + psf = np.full(psf_shape, np.nan) + else: + center = cell_grid.bbox_of(cell_ij).to_legacy().getCenter() + psf = exposure_psf.computeKernelImage(center).array + if psf.shape != psf_shape: + raise ValueError(f"Exposure PSF image at {center=} has shape={psf.shape} != {psf_shape=}") + row_psfs.append(psf) + cell_psfs.append(row_psfs) + + cell_psfs = CellPointSpreadFunction(bounds=cell_grid_bounds, array=np.array(cell_psfs)) + # CellCoadds only take TractFrame for now + sky_projection = SkyProjection.from_legacy( + exposure.wcs, + TractFrame(skymap="dummy", tract=0, bbox=box), + pixel_bounds=box, + ) + coadd_contributions = CoaddProvenance.make_empty_contribution_table(n_rows=1) + coadd_inputs = CoaddProvenance.make_empty_contribution_table(n_rows=1) + patch = PatchDefinition( + id=0, + index=YX(0, 0), + inner_bbox=box, + cells=cell_grid, + ) + + cell_coadd = CellCoadd( + image=Image.from_legacy(exposure.image, unit=image_unit), + variance=Image.from_legacy(exposure.variance), + mask=Mask.from_legacy(exposure.mask), + patch=patch, + provenance=CoaddProvenance(inputs=coadd_inputs, contributions=coadd_contributions), + psf=cell_psfs, + sky_projection=sky_projection, + **kwargs, + ) + return cell_coadd + + def make_test_exposure(): """Make a test exposure with a PSF attached and stars placed randomly. From 2d9149741ee478e2fa0e63b924fd370eb203394d Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Thu, 27 Aug 2026 17:22:33 -0700 Subject: [PATCH 09/23] Add unit test for cell coadd injection --- python/lsst/source/injection/inject_engine.py | 15 +++++++--- tests/test_inject_engine.py | 28 ++++++++++++++++--- 2 files changed, 35 insertions(+), 8 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index f281411..d7c8e6a 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -771,7 +771,7 @@ def inject_galsim_objects_into_exposure( raise TypeError(f"Unsupported {type(exposure)=}") if is_cell: cell_coadd = exposure - psf = exposure.psf + cell_psf = exposure.psf exposure = cell_coadd.to_legacy() is_coadd = True else: @@ -782,6 +782,7 @@ def inject_galsim_objects_into_exposure( exposure.mask.addMaskPlane(mask_plane_core_name) wcs = exposure.getWcs() bbox = exposure.getBBox() + psf = exposure.getPsf() if logger: logger.info( "Adding %s and %s mask planes to the exposure.", @@ -808,10 +809,16 @@ def inject_galsim_objects_into_exposure( cen_bbox = bbox_cell.to_legacy().getCenter() cen_x, cen_y = cen_bbox try: - psf_array = cell_coadd.psf.compute_kernel_image(x=cen_x, y=cen_y) - all_bounds.append((bbox_cell, (cell_ij, cen_x, cen_y))) + psf_array = cell_psf.compute_kernel_image(x=cen_x, y=cen_y).array + bounds_cell = galsim.BoundsI( + bbox_cell.min.x, + bbox_cell.max.x, + bbox_cell.min.y, + bbox_cell.max.y, + ) + all_bounds.append((bounds_cell, (cell_ij, cen_x, cen_y))) if fallback_psf is None: - mat = wcs.linearizePixelToSky(cen_bbox).getMatrix() + mat = wcs.linearizePixelToSky(cen_bbox, arcseconds).getMatrix() galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) fallback_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) except BoundsError: diff --git a/tests/test_inject_engine.py b/tests/test_inject_engine.py index 0627790..836202b 100644 --- a/tests/test_inject_engine.py +++ b/tests/test_inject_engine.py @@ -29,6 +29,8 @@ import lsst.utils.tests from lsst.geom import Point2D, SpherePoint, degrees +from lsst.images import YX +from lsst.images.cells import CellIJ from lsst.source.injection.inject_engine import ( generate_galsim_objects, get_gain_map, @@ -36,7 +38,11 @@ inject_galsim_objects_into_exposure, make_galsim_object, ) -from lsst.source.injection.utils.test_utils import make_test_exposure, make_test_injection_catalog +from lsst.source.injection.utils.test_utils import ( + make_test_cell_coadd, + make_test_exposure, + make_test_injection_catalog, +) from lsst.utils.tests import TestCase @@ -51,6 +57,13 @@ def setUp(self): GalSim objects intended for injection. """ self.exposure = make_test_exposure() + self.cell_coadd = make_test_cell_coadd( + exposure=self.exposure, + cell_shape=YX(x=35, y=35), + psf_shape=(33, 33), + missing={CellIJ(1, 1)}, + band="r", + ) self.injection_catalog = make_test_injection_catalog( self.exposure.getWcs(), self.exposure.getBBox(), @@ -65,6 +78,7 @@ def setUp(self): def tearDown(self): del self.exposure + del self.cell_coadd del self.injection_catalog del self.galsim_objects @@ -137,10 +151,16 @@ def test_get_gain_map_no_valid_pixels(self): self.assertTrue(np.all(np.isfinite(gain_map.array))) self.assertTrue(np.all(gain_map.array > 0)) + def test_inject_galsim_objects_into_cell_coadd(self): + self._test_inject_galsim_objects_into_exposure(self.cell_coadd, True) + def test_inject_galsim_objects_into_exposure(self): - flux0 = np.sum(self.exposure.image.array) + self._test_inject_galsim_objects_into_exposure(self.exposure, False) + + def _test_inject_galsim_objects_into_exposure(self, exposure, is_cell: bool = True): + flux0 = np.sum(exposure.image.array) injected_outputs = inject_galsim_objects_into_exposure( - exposure=self.exposure, + exposure=exposure, objects=self.galsim_objects, mask_plane_name="INJECTED", calib_flux_radius=12.0, @@ -150,7 +170,7 @@ def test_inject_galsim_objects_into_exposure(self): pc = self.exposure.getPhotoCalib() inst_fluxes = [float(pc.magnitudeToInstFlux(mag)) for mag in self.injection_catalog["mag"]] self.assertAlmostEqual( - np.sum(self.exposure.image.array) - flux0, + np.sum(exposure.image.array) - flux0, np.sum(inst_fluxes), delta=0.00015 * np.sum(inst_fluxes), ) From 62e3c9f684d85fddc0e37eafd3b8e61c9f43bc02 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Mon, 31 Aug 2026 16:58:09 -0700 Subject: [PATCH 10/23] Normalize PSF before applying aperture correction --- python/lsst/source/injection/inject_engine.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index d7c8e6a..b4329fa 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -624,6 +624,10 @@ def _get_galsim_psf(*, psf, pixel_coords, bbox, calib_flux_radius, galsim_wcs, s if psf_array is not None: # Compute the aperture corrected PSF interpolated image. + # Normalize first + # TODO: Should also deal with negative pixel values here as they are + # not valid for a PSF. + psf_array /= np.sum(psf_array) # TODO: Review whether we want to at least make this optional aperture_correction = psf.computeApertureFlux(calib_flux_radius, psf.getAveragePosition()) psf_array /= aperture_correction From 8752c2a66f6d2381f6a791bc93b5624e61fb0da5 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Mon, 31 Aug 2026 17:06:21 -0700 Subject: [PATCH 11/23] Make injection_core_size a kwarg --- python/lsst/source/injection/inject_engine.py | 22 ++++++++++++++----- tests/test_inject_engine.py | 9 ++++++++ 2 files changed, 26 insertions(+), 5 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index b4329fa..5ecc257 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -23,6 +23,7 @@ __all__ = ["generate_galsim_objects", "inject_galsim_objects_into_exposure"] +import numbers import os from collections import Counter from collections.abc import Generator @@ -720,6 +721,7 @@ def inject_galsim_objects_into_exposure( add_noise: bool = True, noise_seed: int = 0, bad_mask_names: list[str] | None = None, + injection_core_size: int = 3, logger: Any | None = None, ) -> tuple[list[int], list[galsim.BoundsI], list[bool], list[bool]]: """Inject sources into given exposure using GalSim. @@ -752,6 +754,9 @@ def inject_galsim_objects_into_exposure( Seed for generating the noise of the first injected object. The seed actually used increases by 1 for each subsequent object, to ensure independent noise realizations. + injection_core_size : `int` + Size of the box around each injected source to mark in the injected + core mask. The value must be a positive odd integer. bad_mask_names : `list[str]`, optional List of mask plane names indicating pixels to ignore when fitting flux vs variance in preparation for variance plane modification. If None, @@ -773,6 +778,13 @@ def inject_galsim_objects_into_exposure( is_cell = isinstance(exposure, CellCoadd) if not is_cell and not isinstance(exposure, ExposureF): raise TypeError(f"Unsupported {type(exposure)=}") + if ( + not isinstance(injection_core_size, numbers.Number) + or not (injection_core_size >= 1) + or not (injection_core_size % 2 == 1) + ): + raise ValueError(f"{injection_core_size=} must be a positive odd integer") + if is_cell: cell_coadd = exposure cell_psf = exposure.psf @@ -892,7 +904,7 @@ def inject_galsim_objects_into_exposure( if draw_size == 0: draw_size = conv.getGoodImageSize(galsim_wcs.minLinearScale()) # type: ignore injection_draw_size = int(draw_size) - injection_core_size = 3 + injection_core_size_obj = injection_core_size if draw_size_max > 0 and injection_draw_size > draw_size_max: if logger: logger.warning( @@ -903,15 +915,15 @@ def inject_galsim_objects_into_exposure( ) injection_draw_size = draw_size_max draw_sizes[i] = injection_draw_size - if injection_core_size > injection_draw_size: + if injection_core_size_obj > injection_draw_size: if logger: logger.debug( "Clipping core size for object at %s from %d to %d pixels.", sky_coords, - injection_core_size, + injection_core_size_obj, injection_draw_size, ) - injection_core_size = injection_draw_size + injection_core_size_obj = injection_draw_size sub_bounds = galsim.BoundsI(posi).withBorder(injection_draw_size // 2) # These bounds may not be the same as what would be derived from @@ -949,7 +961,7 @@ def inject_galsim_objects_into_exposure( add_noise=add_noise, noise_seed=noise_seed, is_coadd=is_coadd, - injection_core_size=injection_core_size, + injection_core_size=injection_core_size_obj, mask_plane_name=mask_plane_name, mask_plane_core_name=mask_plane_core_name, logger=logger, diff --git a/tests/test_inject_engine.py b/tests/test_inject_engine.py index 836202b..32903b4 100644 --- a/tests/test_inject_engine.py +++ b/tests/test_inject_engine.py @@ -25,6 +25,7 @@ import galsim import numpy as np +import pytest from galsim import BoundsI, GSObject import lsst.utils.tests @@ -159,6 +160,13 @@ def test_inject_galsim_objects_into_exposure(self): def _test_inject_galsim_objects_into_exposure(self, exposure, is_cell: bool = True): flux0 = np.sum(exposure.image.array) + for injection_core_size in (None, 0, 1.5): + with pytest.raises(ValueError): + inject_galsim_objects_into_exposure( + exposure=exposure, + objects=(), + injection_core_size=injection_core_size, + ) injected_outputs = inject_galsim_objects_into_exposure( exposure=exposure, objects=self.galsim_objects, @@ -166,6 +174,7 @@ def _test_inject_galsim_objects_into_exposure(self, exposure, is_cell: bool = Tr calib_flux_radius=12.0, draw_size_max=1000, add_noise=False, + injection_core_size=5, ) pc = self.exposure.getPhotoCalib() inst_fluxes = [float(pc.magnitudeToInstFlux(mag)) for mag in self.injection_catalog["mag"]] From 0ff73940c64df677cd758bcef54d55098cd26873 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Mon, 31 Aug 2026 17:09:29 -0700 Subject: [PATCH 12/23] Fix remaining cell injection problems --- python/lsst/source/injection/inject_engine.py | 121 ++++++++++++++---- tests/test_inject_engine.py | 49 +++++-- 2 files changed, 136 insertions(+), 34 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index 5ecc257..f5e06c0 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -37,10 +37,12 @@ from astropy.table import Table from galsim import GalSimFFTSizeError +import lsst.afw.detection as afwDetection from lsst.afw.geom import SkyWcs from lsst.afw.image import ExposureF, PhotoCalib from lsst.geom import Box2I, Point2D, Point2I, SpherePoint, arcseconds, degrees -from lsst.images import BoundsError +from lsst.images import BoundsError, MaskPlane +from lsst.images._mask import _guess_legacy_plane_map from lsst.images.cells import CellCoadd, CellIJ from lsst.pex.exceptions import InvalidParameterError, LogicError @@ -604,9 +606,43 @@ def add_noise_to_galsim_image( return image_template, var_template -def _get_galsim_psf(*, psf, pixel_coords, bbox, calib_flux_radius, galsim_wcs, sky_coords=None): - """Compute the galsim PSF at given coordinates within a bbox.""" - psf_array = None +def _get_galsim_psf( + *, + psf: afwDetection.Psf, + pixel_coords: Point2D, + bbox: Box2I, + calib_flux_radius: float, + galsim_wcs: galsim.JacobianWCS, + sky_coords: SpherePoint | None = None, + apply_aperture_correction: bool = True, +): + """Return the GalSim PSF at given coordinates within a bbox. + + Parameters + ---------- + psf + The coadd PSF object to evaluate the PSF with. + pixel_coords + The pixel coordinates to get the PSF for. + bbox + A bounding box to search for an alternative PSF for. If pixel_coords + lies outside the box and PSF evaluation fails, it will retry using the + nearest point to pixel_coords within the box. + calib_flux_radius + The radius to use to compute the aperture flux. + galsim_wcs + The GalSim WCS to pass to the `galsim.InterpolatedImage` constructor. + sky_coords + The sky coordinates of pixel_coords. Only used for + apply_aperture_correction + If True, the PSF values will be divided by the aperture flux at the + average position of psf. + + Returns + ------- + galsim_psf + The GalSim PSF as a `galsim.InterpolatedImage`. + """ error_args = [] try: psf_array = psf.computeKernelImage(pixel_coords).array @@ -622,19 +658,20 @@ def _get_galsim_psf(*, psf, pixel_coords, bbox, calib_flux_radius, galsim_wcs, s psf_array = psf.computeKernelImage(contained_point).array except InvalidParameterError: error_args = ["Cannot compute PSF for object at %s; flagging and skipping.", sky_coords] - - if psf_array is not None: - # Compute the aperture corrected PSF interpolated image. - # Normalize first - # TODO: Should also deal with negative pixel values here as they are - # not valid for a PSF. - psf_array /= np.sum(psf_array) + return None, error_args + + # Compute the aperture corrected PSF interpolated image. + # Normalize first + # TODO: Should also deal with negative pixel values here as they are + # not valid for a PSF. + psf_array /= np.sum(psf_array) + if apply_aperture_correction: # TODO: Review whether we want to at least make this optional aperture_correction = psf.computeApertureFlux(calib_flux_radius, psf.getAveragePosition()) psf_array /= aperture_correction - psf_array = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) + galsim_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) - return psf_array, error_args + return galsim_psf, error_args def _inject_galsim_object_into_bounds( @@ -787,7 +824,7 @@ def inject_galsim_objects_into_exposure( if is_cell: cell_coadd = exposure - cell_psf = exposure.psf + cell_psf = cell_coadd.psf exposure = cell_coadd.to_legacy() is_coadd = True else: @@ -811,7 +848,7 @@ def inject_galsim_objects_into_exposure( pixel_scale = wcs.getPixelScale(bbox.getCenter()).asArcseconds() if is_cell: - all_bounds = [] + all_bounds: list[tuple[galsim.BoundsI, galsim.InterpolatedImage | None]] = [] cell_grid = cell_coadd.grid n_i, n_j = cell_grid.grid_size.as_tuple() @@ -822,23 +859,32 @@ def inject_galsim_objects_into_exposure( for cell_j in range(n_j): cell_ij = CellIJ(cell_i, cell_j) bbox_cell = cell_grid.bbox_of(cell_ij) - cen_bbox = bbox_cell.to_legacy().getCenter() + bbox_cell_legacy = bbox_cell.to_legacy() + cen_bbox = bbox_cell_legacy.getCenter() cen_x, cen_y = cen_bbox try: - psf_array = cell_psf.compute_kernel_image(x=cen_x, y=cen_y).array + if cell_ij in cell_coadd.bounds.missing: + continue + galsim_psf = cell_psf.compute_kernel_image(x=cen_x, y=cen_y).array + mat = wcs.linearizePixelToSky(wcs.pixelToSky(x=cen_x, y=cen_y), arcseconds).getMatrix() + galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) + aperture_correction = psf.computeApertureFlux(calib_flux_radius, cen_bbox) + galsim_psf = galsim.InterpolatedImage( + galsim.Image(galsim_psf / (np.sum(galsim_psf) * aperture_correction)), + wcs=galsim_wcs, + ) bounds_cell = galsim.BoundsI( bbox_cell.min.x, bbox_cell.max.x, bbox_cell.min.y, bbox_cell.max.y, ) - all_bounds.append((bounds_cell, (cell_ij, cen_x, cen_y))) + all_bounds.append((bounds_cell, galsim_psf)) if fallback_psf is None: - mat = wcs.linearizePixelToSky(cen_bbox, arcseconds).getMatrix() - galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) - fallback_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) + fallback_psf = galsim_psf except BoundsError: - # Assume this cell has a bad PSF/data and give up + # TODO: Should this ever happen? If it does, it probably + # means that cell_coadd.bounds.missing is wrong. pass else: all_bounds = [(full_bounds, None)] @@ -937,13 +983,17 @@ def inject_galsim_objects_into_exposure( area_to_inject = 0 area_injected = 0 bounds_fft_size_errors = {} - for idx_bound, (injection_bounds, bounds_info) in enumerate(all_bounds): + for idx_bound, (injection_bounds, galsim_psf_bounds) in enumerate(all_bounds): object_common_bounds = injection_bounds & sub_bounds # Inject the source if there is any overlap. if (common_area := object_common_bounds.area()) > 0: area_to_inject += common_area - if not is_cell: + if is_cell: + # Use the per-cell PSF + assert galsim_psf_bounds is not None + conv = galsim.Convolve(object, galsim_psf_bounds) + else: common_bounds[i] = object_common_bounds # type: ignore try: _inject_galsim_object_into_bounds( @@ -982,4 +1032,27 @@ def inject_galsim_objects_into_exposure( bounds_fft_size_errors, ) + # lsst.images and afw masks are incompatible due to dtypes alone, and so + # they have to be converted explicitly (also to add plane descriptions) + if is_cell: + plane_map = _guess_legacy_plane_map(exposure.mask.getMaskPlaneDict()) + plane_map.update( + { + mask_plane_name: MaskPlane( + name=mask_plane_name, + description="Pixel had at least one source injected over it", + ), + mask_plane_core_name: MaskPlane( + name=mask_plane_core_name, + description=f"Pixel center is within the core region" + f" ({injection_core_size}x{injection_core_size} box) of an injected source", + ), + } + ) + converted = cell_coadd.mask.from_legacy( + exposure.mask, + plane_map=plane_map, + sky_projection=cell_coadd.sky_projection, + ) + cell_coadd.mask = converted return draw_sizes, common_bounds, fft_size_errors, psf_compute_errors diff --git a/tests/test_inject_engine.py b/tests/test_inject_engine.py index 32903b4..a07b390 100644 --- a/tests/test_inject_engine.py +++ b/tests/test_inject_engine.py @@ -58,17 +58,33 @@ def setUp(self): GalSim objects intended for injection. """ self.exposure = make_test_exposure() + # Mark a cell as missing to test that it's ignored, although the data + # in it may still actually be fine + self.cell_bad = CellIJ(1, 1) self.cell_coadd = make_test_cell_coadd( exposure=self.exposure, cell_shape=YX(x=35, y=35), psf_shape=(33, 33), - missing={CellIJ(1, 1)}, + missing={self.cell_bad}, band="r", ) self.injection_catalog = make_test_injection_catalog( self.exposure.getWcs(), self.exposure.getBBox(), ) + cen_cell_bad = self.exposure.wcs.pixelToSky( + self.cell_coadd.grid.bbox_of(CellIJ(1, 1)).to_legacy().getCenter() + ) + row_last = self.injection_catalog[-1] + self.injection_catalog.add_row( + { + "ra": cen_cell_bad.getRa().asDegrees(), + "dec": cen_cell_bad.getDec().asDegrees(), + "mag": row_last["mag"], + "source_type": row_last["source_type"], + } + ) + self.galsim_objects = generate_galsim_objects( injection_catalog=self.injection_catalog, photo_calib=self.exposure.photoCalib, @@ -76,12 +92,18 @@ def setUp(self): fits_alignment="wcs", stamp_prefix="", ) + self.photoCalib = self.exposure.getPhotoCalib() + self.inst_fluxes = [ + float(self.photoCalib.magnitudeToInstFlux(mag)) for mag in self.injection_catalog["mag"] + ] def tearDown(self): del self.exposure del self.cell_coadd del self.injection_catalog del self.galsim_objects + del self.photoCalib + del self.inst_fluxes def test_make_galsim_object(self): source_data = self.injection_catalog[0] @@ -176,19 +198,26 @@ def _test_inject_galsim_objects_into_exposure(self, exposure, is_cell: bool = Tr add_noise=False, injection_core_size=5, ) - pc = self.exposure.getPhotoCalib() - inst_fluxes = [float(pc.magnitudeToInstFlux(mag)) for mag in self.injection_catalog["mag"]] + draw_sizes, common_bounds, fft_size_errors, psf_compute_errors = injected_outputs self.assertAlmostEqual( np.sum(exposure.image.array) - flux0, - np.sum(inst_fluxes), - delta=0.00015 * np.sum(inst_fluxes), + np.sum(self.inst_fluxes), + delta=0.00015 * np.sum(self.inst_fluxes), ) - self.assertEqual(len(injected_outputs[0]), len(self.injection_catalog["ra"])) + self.assertEqual(len(draw_sizes), len(self.injection_catalog["ra"])) self.assertTrue(all(isinstance(injected_output, list) for injected_output in injected_outputs)) - self.assertTrue(all(isinstance(item, int) for item in injected_outputs[0])) # draw sizes - self.assertTrue(all(isinstance(item, BoundsI) for item in injected_outputs[1])) # common bounds - self.assertTrue(all(isinstance(item, bool) for item in injected_outputs[2])) # FFT size errors - self.assertTrue(all(isinstance(item, bool) for item in injected_outputs[3])) # PSF compute errors + self.assertTrue(all(isinstance(item, int) for item in draw_sizes)) + self.assertTrue(all(isinstance(item, BoundsI) for item in common_bounds)) # common bounds + self.assertTrue(all(isinstance(item, bool) for item in fft_size_errors)) # FFT size errors + self.assertTrue(all(isinstance(item, bool) for item in psf_compute_errors)) # PSF compute errors + mask_dict = exposure.mask.schema if is_cell else exposure.mask.getMaskPlaneDict() + assert "INJECTED" in mask_dict + assert "INJECTED_CORE" in mask_dict + # Non-cell coadds shouldn't fail to inject here + # Cell coadds actually should skip this object, but as it is, + # the draw_size is the intended width/length of the injection box, + # not the number of pixels that were actually injected into. + assert draw_sizes[-1] > 0 class MemoryTestCase(lsst.utils.tests.MemoryTestCase): From 4e75b1e0ae99f98061ad4b894757389ab0cb7570 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Tue, 15 Sep 2026 13:56:48 -0700 Subject: [PATCH 13/23] Return empty shear dict if values are all masekd --- python/lsst/source/injection/inject_engine.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index f5e06c0..9b90498 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -131,12 +131,17 @@ def get_shear_data( Dictionary of source data to pass to the Shear constructor. """ shear_params = set(shear_attributes) & set(source_data.keys()) - shear_data = {} - for shear_param in shear_params: - if shear_param == "beta": - shear_data.update({shear_param: source_data[shear_param] * galsim.degrees}) - else: - shear_data.update({shear_param: source_data[shear_param]}) + shear_data = {shear_param: source_data[shear_param] for shear_param in shear_params} + # Return an empty dict if all of the values are masked + # This should only be the case for DeltaFunction, but validation of the + # return is left for the caller of this function. + # While None might perhaps have been a better choice here, changing the + # return type is unnecessary given that galsim.Shear can be initialized + # with no args and correctly returns a trivial shear. + if all(np.ma.is_masked(x) for x in shear_data.values()): + return {} + if (beta := shear_data.get("beta")) is not None: + shear_data["beta"] = beta * galsim.degrees return shear_data From a3a10213b3f94783392f028076c606f1325e6d64 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Tue, 15 Sep 2026 14:00:06 -0700 Subject: [PATCH 14/23] Implement cell injection using same PSF for identical-input cells --- python/lsst/source/injection/inject_engine.py | 79 +++++++++++++++---- 1 file changed, 63 insertions(+), 16 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index 9b90498..01a9be5 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -36,6 +36,8 @@ from astropy.io import fits from astropy.table import Table from galsim import GalSimFFTSizeError +from pydantic import ConfigDict +from pydantic.dataclasses import dataclass import lsst.afw.detection as afwDetection from lsst.afw.geom import SkyWcs @@ -679,6 +681,13 @@ def _get_galsim_psf( return galsim_psf, error_args +@dataclass(config=ConfigDict(arbitrary_types_allowed=True), kw_only=True) +class CellBoundsInfo: + aperture_correction: float = 1.0 + cell_ij: CellIJ + galsim_psf: galsim.InterpolatedImage + + def _inject_galsim_object_into_bounds( *, convolved_object, @@ -853,11 +862,13 @@ def inject_galsim_objects_into_exposure( pixel_scale = wcs.getPixelScale(bbox.getCenter()).asArcseconds() if is_cell: - all_bounds: list[tuple[galsim.BoundsI, galsim.InterpolatedImage | None]] = [] + all_bounds: list[tuple[galsim.BoundsI, CellBoundsInfo | None]] = [] cell_grid = cell_coadd.grid + contributions = cell_coadd.provenance.contributions n_i, n_j = cell_grid.grid_size.as_tuple() fallback_psf = None + cell_inputs: dict[CellIJ, str] = {} # Iterate over all cells and find the ones with good PSFs for cell_i in range(n_i): @@ -870,6 +881,15 @@ def inject_galsim_objects_into_exposure( try: if cell_ij in cell_coadd.bounds.missing: continue + is_cell_input = (contributions["cell_i"] == cell_i) & (contributions["cell_j"] == cell_j) + # Return a string of sorted visits to efficiently count + # unique combinations + cell_inputs[cell_ij] = ",".join( + str(visit) for visit in sorted(set(contributions["visit"][is_cell_input])) + ) + + # Use the native cell coadd PSF methods rather than + # _get_galsim_psf which uses the to_legacy PSF galsim_psf = cell_psf.compute_kernel_image(x=cen_x, y=cen_y).array mat = wcs.linearizePixelToSky(wcs.pixelToSky(x=cen_x, y=cen_y), arcseconds).getMatrix() galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) @@ -884,13 +904,24 @@ def inject_galsim_objects_into_exposure( bbox_cell.min.y, bbox_cell.max.y, ) - all_bounds.append((bounds_cell, galsim_psf)) + bounds_info = CellBoundsInfo( + aperture_correction=aperture_correction, + cell_ij=cell_ij, + galsim_psf=galsim_psf, + ) + all_bounds.append((bounds_cell, bounds_info)) if fallback_psf is None: fallback_psf = galsim_psf except BoundsError: # TODO: Should this ever happen? If it does, it probably # means that cell_coadd.bounds.missing is wrong. pass + cell_visits, cell_visits_index = np.unique(list(cell_inputs.values()), return_inverse=True) + # Store the index of the cell_visits to efficiently check if + # they're identical + cell_inputs_indices = dict(zip(cell_inputs.keys(), cell_visits_index, strict=True)) + # Keep a dict for PSF/bounds retrieval prior to iteration over cells + cell_bounds = {bounds_info.cell_ij: bounds_info for _, bounds_info in all_bounds} else: all_bounds = [(full_bounds, None)] @@ -922,15 +953,6 @@ def inject_galsim_objects_into_exposure( if galsim_pixel_scale < pixel_scale / 2 or galsim_pixel_scale > pixel_scale * 2: continue - # Get the PSF at the centroid of the object - galsim_psf_centroid, psf_error_args = _get_galsim_psf( - psf=psf, - pixel_coords=pixel_coords, - bbox=bbox, - calib_flux_radius=calib_flux_radius, - galsim_wcs=galsim_wcs, - sky_coords=sky_coords, - ) # Get the fallback PSF (if available) to compute the injection box. # Even if the PSF at the centroid of an object is not valid, nearby # cells might be fine. For example, cells with saturated stars can @@ -941,10 +963,26 @@ def inject_galsim_objects_into_exposure( # land outside the exposure's bbox. Non-cell coadds and visits could # also do something similar by searching for the nearest valid PSF, # but that's not as trivial as picking an "average" cell. - if galsim_psf_centroid is None: - if is_cell: + if is_cell: + cell_ij_cen = cell_grid.index_of(x=pixel_coords.x, y=pixel_coords.y) + galsim_psf_centroid = cell_bounds.get(cell_ij_cen) + if galsim_psf_centroid is None: galsim_psf_centroid, psf_error_args = fallback_psf, [] + cell_psf_index_cen = None else: + galsim_psf_centroid = galsim_psf_centroid.galsim_psf + cell_psf_index_cen = cell_inputs_indices[cell_ij_cen] + else: + # Get the PSF at the centroid of the object + galsim_psf_centroid, psf_error_args = _get_galsim_psf( + psf=psf, + pixel_coords=pixel_coords, + bbox=bbox, + calib_flux_radius=calib_flux_radius, + galsim_wcs=galsim_wcs, + sky_coords=sky_coords, + ) + if galsim_psf_centroid is None: psf_compute_errors[i] = True if logger: logger.debug(*psf_error_args) @@ -988,7 +1026,7 @@ def inject_galsim_objects_into_exposure( area_to_inject = 0 area_injected = 0 bounds_fft_size_errors = {} - for idx_bound, (injection_bounds, galsim_psf_bounds) in enumerate(all_bounds): + for idx_bound, (injection_bounds, bounds_info) in enumerate(all_bounds): object_common_bounds = injection_bounds & sub_bounds # Inject the source if there is any overlap. @@ -996,8 +1034,17 @@ def inject_galsim_objects_into_exposure( area_to_inject += common_area if is_cell: # Use the per-cell PSF - assert galsim_psf_bounds is not None - conv = galsim.Convolve(object, galsim_psf_bounds) + galsim_psf = bounds_info.galsim_psf + assert galsim_psf is not None + if bounds_info.cell_ij == cell_ij_cen: + # These have to be one and the same if the cell is not + # missing, which it can't be if it's in all_bounds + assert galsim_psf is galsim_psf_centroid + elif cell_psf_index_cen == cell_inputs_indices[bounds_info.cell_ij]: + # Use the PSF from the cell containing the centroid + assert galsim_psf is not galsim_psf_centroid + galsim_psf = galsim_psf_centroid + conv = galsim.Convolve(object, galsim_psf) else: common_bounds[i] = object_common_bounds # type: ignore try: From 78cbcbc82eba9a92af624418fab8a45a63e3a9fe Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Tue, 15 Sep 2026 14:14:50 -0700 Subject: [PATCH 15/23] Add galaxy to injection test --- tests/test_inject_engine.py | 39 +++++++++++++++++++++++++++++++------ 1 file changed, 33 insertions(+), 6 deletions(-) diff --git a/tests/test_inject_engine.py b/tests/test_inject_engine.py index a07b390..4ea2aae 100644 --- a/tests/test_inject_engine.py +++ b/tests/test_inject_engine.py @@ -68,15 +68,15 @@ def setUp(self): missing={self.cell_bad}, band="r", ) - self.injection_catalog = make_test_injection_catalog( - self.exposure.getWcs(), - self.exposure.getBBox(), - ) cen_cell_bad = self.exposure.wcs.pixelToSky( self.cell_coadd.grid.bbox_of(CellIJ(1, 1)).to_legacy().getCenter() ) - row_last = self.injection_catalog[-1] - self.injection_catalog.add_row( + injection_catalog = make_test_injection_catalog( + self.exposure.getWcs(), + self.exposure.getBBox(), + ) + row_last = injection_catalog[-1] + injection_catalog.add_row( { "ra": cen_cell_bad.getRa().asDegrees(), "dec": cen_cell_bad.getDec().asDegrees(), @@ -84,6 +84,33 @@ def setUp(self): "source_type": row_last["source_type"], } ) + # Add one galaxy that should span multiple cells + self.cell_galaxy = CellIJ(1, 1) + cen_cell_galaxy = self.exposure.wcs.pixelToSky( + self.cell_coadd.grid.bbox_of(CellIJ(1, 1)).to_legacy().getCenter() + ) + n_rows = len(injection_catalog) + columns_sersic = ("n", "half_light_radius", "q", "beta") + injection_catalog.add_columns( + cols=tuple( + np.ma.masked_array(data=np.zeros(n_rows, dtype=float), mask=np.ones(n_rows, dtype=bool)) + for _ in range(len(columns_sersic)) + ), + names=columns_sersic, + ) + injection_catalog.add_row( + { + "ra": cen_cell_galaxy.getRa().asDegrees(), + "dec": cen_cell_galaxy.getDec().asDegrees(), + "mag": row_last["mag"], + "source_type": "Sersic", + "n": 4.0, + "half_light_radius": 20.0, + "q": 0.8, + "beta": 15.0, + } + ) + self.injection_catalog = injection_catalog self.galsim_objects = generate_galsim_objects( injection_catalog=self.injection_catalog, From 551a66fd143e8da787c1cf294d419a2b8152e873 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Fri, 18 Sep 2026 14:46:30 -0700 Subject: [PATCH 16/23] Add CellCoadd I/O support to CoaddInjectTask --- python/lsst/source/injection/inject_base.py | 87 +++++++++++++------- python/lsst/source/injection/inject_coadd.py | 27 +++++- 2 files changed, 82 insertions(+), 32 deletions(-) diff --git a/python/lsst/source/injection/inject_base.py b/python/lsst/source/injection/inject_base.py index ca982d1..b40a688 100644 --- a/python/lsst/source/injection/inject_base.py +++ b/python/lsst/source/injection/inject_base.py @@ -32,8 +32,13 @@ from astropy.table import Table, hstack, vstack from astropy.units import Quantity, UnitConversionError +from lsst.afw.detection import Psf +from lsst.afw.geom import SkyWcs +from lsst.afw.image import ExposureF, PhotoCalib from lsst.afw.image.exposure.exposureUtils import bbox_contains_sky_coords from lsst.geom import Point2D +from lsst.images.cells import CellCoadd +from lsst.images.fields import BaseField from lsst.pex.config import ChoiceField, Field, ListField from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct from lsst.pipe.base.connectionTypes import PrerequisiteInput @@ -206,7 +211,14 @@ class BaseInjectTask(PipelineTask): _DefaultName = "baseInjectTask" ConfigClass = BaseInjectConfig - def run(self, injection_catalogs, input_exposure, psf, photo_calib, wcs): + def run( + self, + injection_catalogs: list[Table], + input_exposure: ExposureF | CellCoadd, + psf: Psf | None = None, + photo_calib: PhotoCalib | None = None, + wcs: SkyWcs | None = None, + ): """Inject sources into an image. Parameters @@ -214,14 +226,16 @@ def run(self, injection_catalogs, input_exposure, psf, photo_calib, wcs): injection_catalogs : `list` [`astropy.table.Table`] Tract level injection catalogs that potentially cover the named input exposure. - input_exposure : `lsst.afw.image.ExposureF` + input_exposure : `lsst.afw.image.ExposureF` or `lsst.images.CellCoadd` The exposure sources will be injected into. psf: `lsst.meas.algorithms.ImagePsf` - PSF model. + PSF model. Only required if input_exposure is an ExposureF. photo_calib : `lsst.afw.image.PhotoCalib` Photometric calibration used to calibrate injected sources. + Only required if input_exposure is an ExposureF. wcs : `lsst.afw.geom.SkyWcs` WCS used to calibrate injected sources. + Only required if input_exposure is an ExposureF. Returns ------- @@ -231,14 +245,24 @@ def run(self, injection_catalogs, input_exposure, psf, photo_calib, wcs): """ self.config = cast(BaseInjectConfig, self.config) - # Attach potential externally calibrated datasets to input_exposure. - # Keep originals so we can reset at the end. - original_psf = input_exposure.getPsf() - original_photo_calib = input_exposure.getPhotoCalib() - original_wcs = input_exposure.getWcs() - input_exposure.setPsf(psf) - input_exposure.setPhotoCalib(photo_calib) - input_exposure.setWcs(wcs) + if is_exposure := isinstance(input_exposure, ExposureF): + # Attach potential externally calibrated datasets to + # input_exposure. Keep originals to reset at the end. + original_psf = input_exposure.getPsf() + original_photo_calib = input_exposure.getPhotoCalib() + original_wcs = input_exposure.getWcs() + input_exposure.setPsf(psf) + input_exposure.setPhotoCalib(photo_calib) + input_exposure.setWcs(wcs) + bbox = input_exposure.getBBox() + photo_calib = input_exposure.getPhotoCalib() + elif isinstance(input_exposure, CellCoadd): + if wcs is None: + wcs = input_exposure.sky_projection.to_legacy() + bbox = input_exposure.bbox.to_legacy() + photo_calib = BaseField.make_legacy_photo_calib(input_exposure.unit) + else: + raise ValueError(f"Unsupported {type(input_exposure)=} for injection") # Make empty table if none supplied to support process_all_data_ids. if len(injection_catalogs) == 0: @@ -274,11 +298,11 @@ def run(self, injection_catalogs, input_exposure, psf, photo_calib, wcs): injection_catalog = self._standardize_columns( injection_catalog, column_mapping, - input_exposure.getWcs().getPixelScale(input_exposure.getBBox().getCenter()).asArcseconds(), + wcs.getPixelScale(bbox.getCenter()).asArcseconds(), ) # Clean the injection catalog of sources which are not injectable. - injection_catalog = self._clean_sources(injection_catalog, input_exposure) + injection_catalog = self._clean_sources(injection_catalog, input_exposure, bbox=bbox, wcs=wcs) # Injection binary flag lookup dictionary. binary_flags = { @@ -379,17 +403,20 @@ def run(self, injection_catalogs, input_exposure, psf, photo_calib, wcs): ) # Restore original input_exposure calibrated data. - input_exposure.setPsf(original_psf) - input_exposure.setPhotoCalib(original_photo_calib) - input_exposure.setWcs(original_wcs) - - # Add injection provenance and injection flags metadata. - metadata = input_exposure.getMetadata() - input_dataset_type = self.config.connections.input_exposure.format(**self.config.connections.toDict()) - metadata.set("INJECTED", input_dataset_type, "Initial source injection dataset type") - input_exposure.getInfo().setVisitInfo(input_exposure.visitInfo.copyWith(hasSimulatedContent=True)) - for flag, value in sorted(binary_flags.items(), key=lambda item: item[1]): - injection_catalog.meta[flag] = value + if is_exposure: + input_exposure.setPsf(original_psf) + input_exposure.setPhotoCalib(original_photo_calib) + input_exposure.setWcs(original_wcs) + + # Add injection provenance and injection flags metadata. + metadata = input_exposure.getMetadata() + input_dataset_type = self.config.connections.input_exposure.format( + **self.config.connections.toDict() + ) + metadata.set("INJECTED", input_dataset_type, "Initial source injection dataset type") + input_exposure.getInfo().setVisitInfo(input_exposure.visitInfo.copyWith(hasSimulatedContent=True)) + for flag, value in sorted(binary_flags.items(), key=lambda item: item[1]): + injection_catalog.meta[flag] = value output_struct = Struct(output_exposure=input_exposure, output_catalog=injection_catalog) return output_struct @@ -494,7 +521,7 @@ def _standardize_columns(self, injection_catalog, column_mapping, pixel_scale): pass return Table(injection_catalog) - def _clean_sources(self, injection_catalog, input_exposure): + def _clean_sources(self, injection_catalog, input_exposure, bbox=None, wcs=None): """Clean the injection catalog of sources which are not injectable. This method will remove sources which are not injectable for a variety @@ -522,6 +549,12 @@ def _clean_sources(self, injection_catalog, input_exposure): """ self.config = cast(BaseInjectConfig, self.config) + if wcs is None: + wcs = input_exposure.getWcs() + + if bbox is None: + bbox = input_exposure.getBBox() + # Exit early if there are no sources to inject. if len(injection_catalog) == 0: self.log.info("Catalog cleaning not applied to empty injection catalog.") @@ -530,13 +563,12 @@ def _clean_sources(self, injection_catalog, input_exposure): sources_to_keep = np.ones(len(injection_catalog), dtype=bool) # Determine centroids and remove sources outside the padded bbox. - wcs = input_exposure.getWcs() has_sky = {"ra", "dec"} <= set(injection_catalog.columns) has_pixel = {"x", "y"} <= set(injection_catalog.columns) # Input catalog must contain either RA/Dec OR x/y. # If only x/y given, RA/Dec will be calculated. if not has_sky and has_pixel: - begin_x, begin_y = input_exposure.getBBox().getBegin() + begin_x, begin_y = bbox.getBegin() ras, decs = wcs.pixelToSkyArray( begin_x + injection_catalog["x"].astype(float), begin_y + injection_catalog["y"].astype(float), @@ -550,7 +582,6 @@ def _clean_sources(self, injection_catalog, input_exposure): elif not has_sky and not has_pixel: self.log.warning("No spatial coordinates found in injection catalog; cannot inject any sources!") if has_sky: - bbox = input_exposure.getBBox() if self.config.trim_padding: bbox.grow(int(self.config.trim_padding)) is_contained = bbox_contains_sky_coords( diff --git a/python/lsst/source/injection/inject_coadd.py b/python/lsst/source/injection/inject_coadd.py index 8147b88..fae24b7 100644 --- a/python/lsst/source/injection/inject_coadd.py +++ b/python/lsst/source/injection/inject_coadd.py @@ -23,6 +23,10 @@ __all__ = ["CoaddInjectConnections", "CoaddInjectConfig", "CoaddInjectTask"] +import dataclasses +from typing import cast + +from lsst.pex.config import Field from lsst.pipe.base.connectionTypes import Input, Output from .inject_base import BaseInjectConfig, BaseInjectConnections, BaseInjectTask @@ -56,6 +60,16 @@ class CoaddInjectConnections( dimensions=("skymap", "tract", "patch", "band"), ) + def __init__(self, *, config: CoaddInjectConfig | None = None): + super().__init__(config=config) + + if not self.config: + return + + if self.config.use_cell_coadds: + self.input_exposure = dataclasses.replace(self.input_exposure, storageClass="CellCoadd") + self.output_exposure = dataclasses.replace(self.output_exposure, storageClass="CellCoadd") + class CoaddInjectConfig( # type: ignore [call-arg] BaseInjectConfig, @@ -63,7 +77,7 @@ class CoaddInjectConfig( # type: ignore [call-arg] ): """Coadd-level configuration for source injection tasks.""" - pass + use_cell_coadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?") class CoaddInjectTask(BaseInjectTask): @@ -75,9 +89,14 @@ class CoaddInjectTask(BaseInjectTask): def runQuantum(self, butler_quantum_context, input_refs, output_refs): inputs = butler_quantum_context.get(input_refs) - inputs["psf"] = inputs["input_exposure"].getPsf() - inputs["photo_calib"] = inputs["input_exposure"].getPhotoCalib() - inputs["wcs"] = inputs["input_exposure"].getWcs() + config = cast(self.ConfigClass, self.config) + exposure = inputs["input_exposure"] + if config.use_cell_coadds: + exposure.apply_background(None) + else: + inputs["psf"] = exposure.getPsf() + inputs["photo_calib"] = exposure.getPhotoCalib() + inputs["wcs"] = exposure.getWcs() input_keys = ["injection_catalogs", "input_exposure", "sky_map", "psf", "photo_calib", "wcs"] outputs = self.run(**{key: value for (key, value) in inputs.items() if key in input_keys}) From 70f9c77c01efafb87ca2843acca548bdf559f493 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Fri, 18 Sep 2026 14:52:27 -0700 Subject: [PATCH 17/23] Inject Iterable of objects instead of Generator --- python/lsst/source/injection/inject_engine.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index 01a9be5..ee631e3 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -26,7 +26,7 @@ import numbers import os from collections import Counter -from collections.abc import Generator +from collections.abc import Generator, Iterable from typing import Any import galsim @@ -764,7 +764,7 @@ def _inject_galsim_object_into_bounds( def inject_galsim_objects_into_exposure( exposure: ExposureF | CellCoadd, - objects: Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None], + objects: Iterable[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None], mask_plane_name: str = "INJECTED", calib_flux_radius: float = 12.0, draw_size_max: int = 1000, @@ -781,7 +781,7 @@ def inject_galsim_objects_into_exposure( ---------- exposure : `lsst.afw.image.ExposureF` or `lsst.images.cells.CellCoadd` The exposure to inject synthetic sources into. - objects : `Generator` [`tuple`, None, None] + objects : `Iterable` [`tuple`, None, None] An iterator of tuples that contains (or generates) locations and object surface brightness profiles to inject. The tuples should contain the following elements: `lsst.geom.SpherePoint`, `lsst.geom.Point2D`, From b6baf9000cbf497cdcb438083ddf6c77fd4e0265 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Mon, 21 Sep 2026 09:19:16 -0700 Subject: [PATCH 18/23] Only add core mask in cell containing centroid --- python/lsst/source/injection/inject_engine.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index ee631e3..5d451fd 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -1048,12 +1048,20 @@ def inject_galsim_objects_into_exposure( else: common_bounds[i] = object_common_bounds # type: ignore try: + if is_cell: + # only the cell containing the centroid should have + # the core mask plane set + injection_core_bounds = full_bounds & object_common_bounds.withBorder( + injection_core_size // 2 + ) + else: + injection_core_bounds = full_bounds _inject_galsim_object_into_bounds( convolved_object=conv, position_d=posd, position_i=posi, object_common_bounds=object_common_bounds, - full_bounds=full_bounds, + full_bounds=injection_core_bounds, galsim_image=galsim_image, galsim_variance=galsim_variance, galsim_wcs=galsim_wcs, @@ -1068,7 +1076,7 @@ def inject_galsim_objects_into_exposure( mask_plane_core_name=mask_plane_core_name, logger=logger, ) - area_injected += 0 + area_injected += common_area except GalSimFFTSizeError as err: bounds_fft_size_errors[idx_bound] = err.size if logger and (area_to_inject == 0): From bf496e166710b3531c5c91d9980b184a5a6c3f2e Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Mon, 21 Sep 2026 09:19:58 -0700 Subject: [PATCH 19/23] Expand injection unit test This adds galaxy injection and more thorough testing for missing cells. --- tests/test_inject_engine.py | 174 ++++++++++++++++++++++++++---------- 1 file changed, 126 insertions(+), 48 deletions(-) diff --git a/tests/test_inject_engine.py b/tests/test_inject_engine.py index 4ea2aae..5daf3f7 100644 --- a/tests/test_inject_engine.py +++ b/tests/test_inject_engine.py @@ -23,6 +23,7 @@ import unittest from types import GeneratorType +import astropy.units as u import galsim import numpy as np import pytest @@ -60,20 +61,39 @@ def setUp(self): self.exposure = make_test_exposure() # Mark a cell as missing to test that it's ignored, although the data # in it may still actually be fine - self.cell_bad = CellIJ(1, 1) - self.cell_coadd = make_test_cell_coadd( - exposure=self.exposure, - cell_shape=YX(x=35, y=35), - psf_shape=(33, 33), - missing={self.cell_bad}, - band="r", + self.cell_bad = CellIJ(3, 3) + self.cell_bad_inverse = { + CellIJ(i, j) + for i in range(90) + for j in range(90) + if not (i == self.cell_bad.i and j == self.cell_bad.j) + } + self.cell_coadd, self.cell_coadd_inverse = ( + make_test_cell_coadd( + exposure=self.exposure.clone(), + cell_shape=YX(x=35, y=35), + psf_shape=(33, 33), + missing=missing, + band="r", + ) + for missing in ({self.cell_bad}, self.cell_bad_inverse) + ) + bbox_exp = self.exposure.getBBox() + bbox_cell_bad = self.cell_coadd.grid.bbox_of(self.cell_bad) + self.slice_x_cell_bad = slice( + bbox_cell_bad.start.x - bbox_exp.beginX, + bbox_cell_bad.stop.x - bbox_exp.beginX, ) - cen_cell_bad = self.exposure.wcs.pixelToSky( - self.cell_coadd.grid.bbox_of(CellIJ(1, 1)).to_legacy().getCenter() + self.slice_y_cell_bad = slice( + bbox_cell_bad.start.y - bbox_exp.beginY, + bbox_cell_bad.stop.y - bbox_exp.beginY, ) + self.bbox_cell_bad = bbox_cell_bad + self.bbox_cell_bad_legacy = bbox_cell_bad.to_legacy() + cen_cell_bad = self.exposure.wcs.pixelToSky(self.bbox_cell_bad_legacy.getCenter()) injection_catalog = make_test_injection_catalog( self.exposure.getWcs(), - self.exposure.getBBox(), + bbox_exp, ) row_last = injection_catalog[-1] injection_catalog.add_row( @@ -84,10 +104,10 @@ def setUp(self): "source_type": row_last["source_type"], } ) - # Add one galaxy that should span multiple cells - self.cell_galaxy = CellIJ(1, 1) + # Add one galaxy centered in the "bad" cell + self.cell_galaxy = self.cell_bad cen_cell_galaxy = self.exposure.wcs.pixelToSky( - self.cell_coadd.grid.bbox_of(CellIJ(1, 1)).to_legacy().getCenter() + self.cell_coadd.grid.bbox_of(self.cell_galaxy).to_legacy().getCenter() ) n_rows = len(injection_catalog) columns_sersic = ("n", "half_light_radius", "q", "beta") @@ -102,14 +122,15 @@ def setUp(self): { "ra": cen_cell_galaxy.getRa().asDegrees(), "dec": cen_cell_galaxy.getDec().asDegrees(), - "mag": row_last["mag"], + "mag": injection_catalog["mag"][0], "source_type": "Sersic", - "n": 4.0, - "half_light_radius": 20.0, + "n": 2.0, + "half_light_radius": 10 * 0.036, "q": 0.8, "beta": 15.0, } ) + injection_catalog["half_light_radius"].unit = u.arcsec self.injection_catalog = injection_catalog self.galsim_objects = generate_galsim_objects( @@ -202,13 +223,24 @@ def test_get_gain_map_no_valid_pixels(self): self.assertTrue(np.all(gain_map.array > 0)) def test_inject_galsim_objects_into_cell_coadd(self): - self._test_inject_galsim_objects_into_exposure(self.cell_coadd, True) + self._test_inject_galsim_objects_into_exposure( + self.cell_coadd, + exposure_inverse=self.cell_coadd_inverse, + is_cell=True, + ) def test_inject_galsim_objects_into_exposure(self): - self._test_inject_galsim_objects_into_exposure(self.exposure, False) + self._test_inject_galsim_objects_into_exposure(self.exposure, is_cell=False) - def _test_inject_galsim_objects_into_exposure(self, exposure, is_cell: bool = True): - flux0 = np.sum(exposure.image.array) + def _test_inject_galsim_objects_into_exposure( + self, + exposure, + exposure_inverse=None, + is_cell: bool = True, + ): + array_orig = exposure.image.array.copy() + flux0 = np.sum(array_orig) + cell_bad_orig = array_orig[self.slice_y_cell_bad, self.slice_x_cell_bad] for injection_core_size in (None, 0, 1.5): with pytest.raises(ValueError): inject_galsim_objects_into_exposure( @@ -216,35 +248,81 @@ def _test_inject_galsim_objects_into_exposure(self, exposure, is_cell: bool = Tr objects=(), injection_core_size=injection_core_size, ) - injected_outputs = inject_galsim_objects_into_exposure( - exposure=exposure, - objects=self.galsim_objects, - mask_plane_name="INJECTED", - calib_flux_radius=12.0, - draw_size_max=1000, - add_noise=False, - injection_core_size=5, - ) - draw_sizes, common_bounds, fft_size_errors, psf_compute_errors = injected_outputs - self.assertAlmostEqual( - np.sum(exposure.image.array) - flux0, - np.sum(self.inst_fluxes), - delta=0.00015 * np.sum(self.inst_fluxes), + galsim_objects = self.galsim_objects + if is_cell: + # At least the last object will be needed later + galsim_objects = list(galsim_objects) + + injected_outputs, injected_outputs_inverse = ( + ( + inject_galsim_objects_into_exposure( + exposure=to_inject, + objects=galsim_objects, + mask_plane_name="INJECTED", + calib_flux_radius=12.0, + draw_size_max=1000, + add_noise=False, + injection_core_size=5, + ) + if to_inject is not None + else None + ) + for to_inject in (exposure, exposure_inverse) ) - self.assertEqual(len(draw_sizes), len(self.injection_catalog["ra"])) - self.assertTrue(all(isinstance(injected_output, list) for injected_output in injected_outputs)) - self.assertTrue(all(isinstance(item, int) for item in draw_sizes)) - self.assertTrue(all(isinstance(item, BoundsI) for item in common_bounds)) # common bounds - self.assertTrue(all(isinstance(item, bool) for item in fft_size_errors)) # FFT size errors - self.assertTrue(all(isinstance(item, bool) for item in psf_compute_errors)) # PSF compute errors - mask_dict = exposure.mask.schema if is_cell else exposure.mask.getMaskPlaneDict() - assert "INJECTED" in mask_dict - assert "INJECTED_CORE" in mask_dict - # Non-cell coadds shouldn't fail to inject here - # Cell coadds actually should skip this object, but as it is, - # the draw_size is the intended width/length of the injection box, - # not the number of pixels that were actually injected into. - assert draw_sizes[-1] > 0 + + cell_bad = exposure.image.array[self.slice_y_cell_bad, self.slice_x_cell_bad] + mask_cell_bad = exposure.mask[self.bbox_cell_bad if is_cell else self.bbox_cell_bad_legacy] + + check_total_flux = True + if is_cell: + # Test that nothing was injected in the bad cell + np.testing.assert_array_equal(cell_bad_orig, cell_bad) + np.testing.assert_array_equal(mask_cell_bad.get("INJECTED"), False) + np.testing.assert_array_equal(mask_cell_bad.get("INJECTED_CORE"), False) + if exposure_inverse is not None: + cell_bad_inverse = exposure_inverse.image.array[self.slice_y_cell_bad, self.slice_x_cell_bad] + mask_cell_bad_inverse = exposure_inverse.mask[self.bbox_cell_bad] + np.testing.assert_array_less(cell_bad_orig, cell_bad_inverse) + np.testing.assert_array_equal(mask_cell_bad_inverse.get("INJECTED"), True) + # Default core size is 3x3 and there should be at least one + # injected core from the galaxy. Others may exist depending on + # the random seed used + self.assertGreaterEqual(np.sum(mask_cell_bad_inverse.get("INJECTED_CORE")), 9) + exposure.image[self.bbox_cell_bad].array = exposure_inverse.image[self.bbox_cell_bad].array + else: + check_total_flux = False + else: + # Test that the entire good cell is filled + # Note this requires that the injected galaxy be large enough, so + # if setUp is modified to inject a smaller object it will fail + np.testing.assert_array_less(cell_bad_orig, cell_bad) + injected_bitmask = mask_cell_bad.getPlaneBitMask("INJECTED") + np.testing.assert_array_equal(mask_cell_bad.array & injected_bitmask, injected_bitmask) + + if check_total_flux: + self.assertAlmostEqual( + np.sum(exposure.image.array) - flux0, + np.sum(self.inst_fluxes), + delta=0.00015 * np.sum(self.inst_fluxes), + ) + + for injected_exposure in (exposure, exposure_inverse): + if injected_exposure is None: + continue + mask_dict = exposure.mask.schema if is_cell else exposure.mask.getMaskPlaneDict() + assert "INJECTED" in mask_dict + assert "INJECTED_CORE" in mask_dict + + for injected_outputs in (injected_outputs, injected_outputs_inverse): + if injected_outputs is None: + continue + draw_sizes, common_bounds, fft_size_errors, psf_compute_errors = injected_outputs + self.assertEqual(len(draw_sizes), len(self.injection_catalog["ra"])) + self.assertTrue(all(isinstance(injected_output, list) for injected_output in injected_outputs)) + self.assertTrue(all(isinstance(item, int) and (item > 0) for item in draw_sizes)) + self.assertTrue(all(isinstance(item, BoundsI) for item in common_bounds)) # common bounds + self.assertTrue(all(item is False for item in fft_size_errors)) # FFT size errors + self.assertTrue(all(isinstance(item, bool) for item in psf_compute_errors)) # PSF compute errors class MemoryTestCase(lsst.utils.tests.MemoryTestCase): From 9360690b3c742895c6369cbebe851fb714de20f9 Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Wed, 7 Oct 2026 20:24:35 -0700 Subject: [PATCH 20/23] Remove gaap per-band flags from matched catalogs These are no longer present in the object table. --- pipelines/match_injected_tract_catalog.yaml | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/pipelines/match_injected_tract_catalog.yaml b/pipelines/match_injected_tract_catalog.yaml index 727f72a..5a68a8d 100644 --- a/pipelines/match_injected_tract_catalog.yaml +++ b/pipelines/match_injected_tract_catalog.yaml @@ -57,11 +57,9 @@ tasks: columns_target_copy.append(f"{band}_gaapFlux_flag") done_gaap = True suffixes = ["Flux", "FluxErr"] - if is_gaap: - suffixes.append("Flux_flag_bigPsf") - elif algo == "cModel": + if algo == "cModel": suffixes.append("_flag") - elif algo != "sersic": + elif (algo != "sersic") and not algo.startswith("gaap"): suffixes.append("Flux_flag") for suffix in suffixes: columns_target_copy.append(f"{band}_{algo}{suffix}") From 515fa90eab6b939fdc751d0f29415993ef4dc6da Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Wed, 7 Oct 2026 22:08:10 -0700 Subject: [PATCH 21/23] Remove gaap3p0 columns from matched catalogs --- pipelines/match_injected_tract_catalog.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pipelines/match_injected_tract_catalog.yaml b/pipelines/match_injected_tract_catalog.yaml index 5a68a8d..7a276be 100644 --- a/pipelines/match_injected_tract_catalog.yaml +++ b/pipelines/match_injected_tract_catalog.yaml @@ -1,7 +1,7 @@ description: Match tract-level injected coadd catalogs to injected object table datasets. parameters: bands_injected: ["u", "g", "r", "i", "z", "y"] - flux_algorithms: ["cModel", "gaap1p0", "gaap3p0", "kron", "psf", "sersic"] + flux_algorithms: ["cModel", "gaap1p0", "kron", "psf", "sersic"] source_type_column: "source_type" tasks: consolidateInjectedCatalogs: From bfb4c7e9e262ae910dbe02f0d2aaaa737755f71c Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Wed, 7 Oct 2026 22:40:52 -0700 Subject: [PATCH 22/23] Remove sersic_x/y columns from matched catalogs --- pipelines/match_injected_tract_catalog.yaml | 1 - 1 file changed, 1 deletion(-) diff --git a/pipelines/match_injected_tract_catalog.yaml b/pipelines/match_injected_tract_catalog.yaml index 7a276be..956f994 100644 --- a/pipelines/match_injected_tract_catalog.yaml +++ b/pipelines/match_injected_tract_catalog.yaml @@ -65,7 +65,6 @@ tasks: columns_target_copy.append(f"{band}_{algo}{suffix}") if "sersic" in parameters.flux_algorithms: columns_target_copy.extend([ - "sersic_x", "sersic_y", "sersic_xErr", "sersic_yErr", "sersic_ra", "sersic_dec", "sersic_raErr", "sersic_decErr", "sersic_index", "sersic_indexErr", "sersic_rho", "sersic_rhoErr", "sersic_reff_x", "sersic_reff_xErr", "sersic_reff_y", "sersic_reff_yErr", From 99410ece71dbbe3cf349662ee744b1f0591d60ce Mon Sep 17 00:00:00 2001 From: Dan Taranu Date: Wed, 7 Oct 2026 22:42:20 -0700 Subject: [PATCH 23/23] Add useful sersic shape columns to matched catalogs --- pipelines/match_injected_tract_catalog.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/pipelines/match_injected_tract_catalog.yaml b/pipelines/match_injected_tract_catalog.yaml index 956f994..b3f8b7c 100644 --- a/pipelines/match_injected_tract_catalog.yaml +++ b/pipelines/match_injected_tract_catalog.yaml @@ -68,6 +68,7 @@ tasks: "sersic_ra", "sersic_dec", "sersic_raErr", "sersic_decErr", "sersic_index", "sersic_indexErr", "sersic_rho", "sersic_rhoErr", "sersic_reff_x", "sersic_reff_xErr", "sersic_reff_y", "sersic_reff_yErr", + "sersic_reff_major", "sersic_reff_minor", "sersic_theta", "sersic_unknown_flag", "sersic_no_data_flag", ]) # TODO: Remove as part of DM-44139