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: diff --git a/pipelines/match_injected_tract_catalog.yaml b/pipelines/match_injected_tract_catalog.yaml index 50953fb..b3f8b7c 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", "kron", "psf", "sersic"] + source_type_column: "source_type" tasks: consolidateInjectedCatalogs: class: lsst.source.injection.utils.ConsolidateInjectedCatalogsTask @@ -19,40 +23,53 @@ 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_" - # 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 + 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 algo == "cModel": + suffixes.append("_flag") + elif (algo != "sersic") and not algo.startswith("gaap"): + 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_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 + config.diff_matched_catalog.columns_ref_mag_to_nJy = columns_ref_mag_to_nJy 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}) diff --git a/python/lsst/source/injection/inject_engine.py b/python/lsst/source/injection/inject_engine.py index 63e8a45..5d451fd 100644 --- a/python/lsst/source/injection/inject_engine.py +++ b/python/lsst/source/injection/inject_engine.py @@ -23,23 +23,34 @@ __all__ = ["generate_galsim_objects", "inject_galsim_objects_into_exposure"] +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 +import galsim.errors as galsim_errors import numpy as np import numpy.ma as ma 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 from lsst.afw.image import ExposureF, PhotoCalib from lsst.geom import Box2I, Point2D, Point2I, SpherePoint, arcseconds, degrees +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 +# 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 @@ -122,12 +133,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 @@ -597,9 +613,158 @@ def add_noise_to_galsim_image( return image_template, var_template +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 + 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] + 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 + galsim_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs) + + 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, + 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 + + 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], + exposure: ExposureF | CellCoadd, + 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, @@ -607,15 +772,16 @@ 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. 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] + 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`, @@ -639,6 +805,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, @@ -657,25 +826,107 @@ 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 ( + 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 = cell_coadd.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() + psf = exposure.getPsf() 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: 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): + for cell_j in range(n_j): + cell_ij = CellIJ(cell_i, cell_j) + bbox_cell = cell_grid.bbox_of(cell_ij) + bbox_cell_legacy = bbox_cell.to_legacy() + cen_bbox = bbox_cell_legacy.getCenter() + cen_x, cen_y = cen_bbox + 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]) + 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, + ) + 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)] + + # 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] = [] @@ -702,39 +953,47 @@ 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) + # 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 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 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: + if galsim_psf_centroid is None: psf_compute_errors[i] = True if logger: - logger.debug("Cannot compute PSF for object at %s; flagging and skipping.", sky_coords) + logger.debug(*psf_error_args) 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) - # 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) - 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( @@ -745,88 +1004,115 @@ 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) - 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: - 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" - ) - 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, - ) - 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: + # 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 is_cell: + # Use the per-cell PSF + 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: + 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=injection_core_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_obj, + mask_plane_name=mask_plane_name, + mask_plane_core_name=mask_plane_core_name, + logger=logger, + ) + area_injected += common_area + 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, + ) + # 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/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 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. diff --git a/tests/test_inject_engine.py b/tests/test_inject_engine.py index 0627790..5daf3f7 100644 --- a/tests/test_inject_engine.py +++ b/tests/test_inject_engine.py @@ -23,12 +23,16 @@ import unittest from types import GeneratorType +import astropy.units as u import galsim import numpy as np +import pytest from galsim import BoundsI, GSObject 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 +40,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,10 +59,80 @@ def setUp(self): GalSim objects intended for injection. """ self.exposure = make_test_exposure() - self.injection_catalog = make_test_injection_catalog( + # 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(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, + ) + 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( + { + "ra": cen_cell_bad.getRa().asDegrees(), + "dec": cen_cell_bad.getDec().asDegrees(), + "mag": row_last["mag"], + "source_type": row_last["source_type"], + } ) + # 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(self.cell_galaxy).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": injection_catalog["mag"][0], + "source_type": "Sersic", + "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( injection_catalog=self.injection_catalog, photo_calib=self.exposure.photoCalib, @@ -62,11 +140,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] @@ -137,29 +222,107 @@ 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_exposure(self): - flux0 = np.sum(self.exposure.image.array) - injected_outputs = inject_galsim_objects_into_exposure( - exposure=self.exposure, - objects=self.galsim_objects, - mask_plane_name="INJECTED", - calib_flux_radius=12.0, - draw_size_max=1000, - add_noise=False, + def test_inject_galsim_objects_into_cell_coadd(self): + self._test_inject_galsim_objects_into_exposure( + self.cell_coadd, + exposure_inverse=self.cell_coadd_inverse, + is_cell=True, ) - 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(inst_fluxes), - delta=0.00015 * np.sum(inst_fluxes), + + def test_inject_galsim_objects_into_exposure(self): + self._test_inject_galsim_objects_into_exposure(self.exposure, is_cell=False) + + 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( + exposure=exposure, + objects=(), + injection_core_size=injection_core_size, + ) + 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(injected_outputs[0]), 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 + + 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):