Coverage for python/lsst/source/injection/inject_base.py: 17%
232 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-11 09:07 +0000
1# This file is part of source_injection.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22from __future__ import annotations
24__all__ = ["_ALLOWED_SOURCE_TYPES", "BaseInjectConnections", "BaseInjectConfig", "BaseInjectTask"]
26from typing import cast
28import galsim
29import numpy as np
30import numpy.ma as ma
31from astropy import units
32from astropy.table import Table, hstack, vstack
33from astropy.units import Quantity, UnitConversionError
35from lsst.afw.detection import Psf
36from lsst.afw.geom import SkyWcs
37from lsst.afw.image import ExposureF, PhotoCalib
38from lsst.afw.image.exposure.exposureUtils import bbox_contains_sky_coords
39from lsst.geom import Point2D
40from lsst.images.cells import CellCoadd
41from lsst.images.fields import BaseField
42from lsst.pex.config import ChoiceField, Field, ListField
43from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
44from lsst.pipe.base.connectionTypes import PrerequisiteInput
46from .inject_engine import generate_galsim_objects, inject_galsim_objects_into_exposure
48_ALLOWED_SOURCE_TYPES = [
49 "Gaussian",
50 "Box",
51 "TopHat",
52 "DeltaFunction",
53 "Airy",
54 "Moffat",
55 "Kolmogorov",
56 "VonKarman",
57 "Exponential",
58 "DeVaucouleurs",
59 "Sersic",
60 "InclinedExponential",
61 "InclinedSersic",
62 "Spergel",
63 "RandomKnots",
64 "Star",
65 "Trail",
66 "Stamp",
67]
70class BaseInjectConnections(
71 PipelineTaskConnections,
72 dimensions=("instrument",),
73 defaultTemplates={
74 "injection_prefix": "injection_",
75 "injected_prefix": "injected_",
76 },
77):
78 """Base connections for source injection tasks."""
80 injection_catalogs = PrerequisiteInput(
81 doc="Set of catalogs of sources to draw inputs from.",
82 name="{injection_prefix}catalog",
83 dimensions=("htm7", "band"),
84 storageClass="ArrowAstropy",
85 minimum=0,
86 multiple=True,
87 )
90class BaseInjectConfig(PipelineTaskConfig, pipelineConnections=BaseInjectConnections):
91 """Base configuration for source injection tasks."""
93 # Catalog manipulation options.
94 process_all_data_ids = Field[bool](
95 doc="If True, all input data IDs will be processed, even those where no synthetic sources were "
96 "identified for injection. In such an eventuality this returns a clone of the input image, renamed "
97 "to the *output_exposure* connection name and with an empty *mask_plane_name* mask plane attached.",
98 default=False,
99 )
100 trim_padding = Field[int](
101 doc="Size of the pixel padding surrounding the image. Only those synthetic sources with a centroid "
102 "falling within the ``image + trim_padding`` region will be considered for source injection.",
103 default=100,
104 optional=True,
105 )
106 selection = Field[str](
107 doc="A string that can be evaluated as a boolean expression to select rows in the input injection "
108 "catalog. To make use of this configuration option, the internal object name ``injection_catalog`` "
109 "must be used. For example, to select all sources with a magnitude in the range 20.0 < mag < 25.0, "
110 "set ``selection=\"(injection_catalog['mag'] > 20.0) & (injection_catalog['mag'] < 25.0)\"``. "
111 "The ``{visit}`` field will be substituted for the current visit ID of the exposure being processed. "
112 "For example, to select only visits that match a user-supplied visit column in the input injection "
113 "catalog, set ``selection=\"np.isin(injection_catalog['visit'], {visit})\"``.",
114 optional=True,
115 )
117 # General configuration options.
118 mask_plane_name = Field[str](
119 doc="Name assigned to the injected mask plane which is attached to the output exposure.",
120 default="INJECTED",
121 )
122 calib_flux_radius = Field[float](
123 doc="Aperture radius (in pixels) that was used to define the calibration for this image+catalog. "
124 "This will be used to produce the correct instrumental fluxes within the radius. "
125 "This value should match that of the field defined in ``slot_CalibFlux_instFlux``.",
126 default=12.0,
127 )
128 fits_alignment = ChoiceField[str]( # type: ignore
129 doc="How should injections from FITS files be aligned?",
130 dtype=str,
131 allowed={
132 "wcs": (
133 "Input image will be transformed such that the local WCS in the FITS header matches the "
134 "local WCS in the target image. I.e., North, East, and angular distances in the input image "
135 "will match North, East, and angular distances in the target image."
136 ),
137 "pixel": (
138 "Input image will **not** be transformed. Up, right, and pixel distances in the input image "
139 "will match up, right and pixel distances in the target image."
140 ),
141 },
142 default="pixel",
143 )
144 stamp_prefix = Field[str](
145 doc="String to prefix to the entries in the *col_stamp* column, for example, a directory path.",
146 default="",
147 )
148 inject_variance = Field[bool](
149 doc="Whether, when injecting flux into the image plane, to inject a corresponding amount of variance "
150 "into the variance plane.",
151 default=True,
152 )
153 add_noise = Field[bool](
154 doc="Whether to randomly vary the injected flux in each pixel by an amount consistent with "
155 "the injected variance.",
156 default=True,
157 )
158 noise_seed = Field[int](
159 doc="Initial seed for random noise generation. This value increments by 1 for each injected "
160 "object, so each object has an independent noise realization.",
161 default=0,
162 )
163 bad_mask_names = ListField[str](
164 doc="List of mask plane names indicating pixels to ignore when fitting flux vs variance in "
165 "preparation for variance plane modification.",
166 default=["BAD", "CR", "CROSSTALK", "INTRP", "NO_DATA", "SAT", "SUSPECT", "UNMASKEDNAN"],
167 )
169 # Custom column names.
170 col_ra = Field[str](
171 doc="Column name for right ascension (in degrees).",
172 default="ra",
173 )
174 col_dec = Field[str](
175 doc="Column name for declination (in degrees).",
176 default="dec",
177 )
178 col_source_type = Field[str](
179 doc="Column name for the source type used in the input catalog. Must match one of the surface "
180 "brightness profiles defined by GalSim.",
181 default="source_type",
182 )
183 col_mag = Field[str](
184 doc="Column name for magnitude.",
185 default="mag",
186 )
187 col_stamp = Field[str](
188 doc="Column name to identify FITS file postage stamps for direct injection. The strings in this "
189 "column will be prefixed with a string given in *stamp_prefix*, to assist in providing the full "
190 "path to a FITS file.",
191 default="stamp",
192 )
193 col_draw_size = Field[str](
194 doc="Column name providing pixel size of the region into which the source profile will be drawn. If "
195 "this column is not provided as an input, the GalSim method ``getGoodImageSize`` will be used "
196 "instead.",
197 default="draw_size",
198 )
199 col_trail_length = Field[str](
200 doc="Column name for specifying a satellite trail length (in pixels).",
201 default="trail_length",
202 )
204 def setDefaults(self):
205 super().setDefaults()
208class BaseInjectTask(PipelineTask):
209 """Base class for injecting sources into images."""
211 _DefaultName = "baseInjectTask"
212 ConfigClass = BaseInjectConfig
214 def run(
215 self,
216 injection_catalogs: list[Table],
217 input_exposure: ExposureF | CellCoadd,
218 psf: Psf | None = None,
219 photo_calib: PhotoCalib | None = None,
220 wcs: SkyWcs | None = None,
221 ):
222 """Inject sources into an image.
224 Parameters
225 ----------
226 injection_catalogs : `list` [`astropy.table.Table`]
227 Tract level injection catalogs that potentially cover the named
228 input exposure.
229 input_exposure : `lsst.afw.image.ExposureF` or `lsst.images.CellCoadd`
230 The exposure sources will be injected into.
231 psf: `lsst.meas.algorithms.ImagePsf`
232 PSF model. Only required if input_exposure is an ExposureF.
233 photo_calib : `lsst.afw.image.PhotoCalib`
234 Photometric calibration used to calibrate injected sources.
235 Only required if input_exposure is an ExposureF.
236 wcs : `lsst.afw.geom.SkyWcs`
237 WCS used to calibrate injected sources.
238 Only required if input_exposure is an ExposureF.
240 Returns
241 -------
242 output_struct : `lsst.pipe.base.Struct`
243 contains : output_exposure : `lsst.afw.image.ExposureF`
244 output_catalog : `lsst.afw.table.SourceCatalog`
245 """
246 self.config = cast(BaseInjectConfig, self.config)
248 if is_exposure := isinstance(input_exposure, ExposureF):
249 # Attach potential externally calibrated datasets to
250 # input_exposure. Keep originals to reset at the end.
251 original_psf = input_exposure.getPsf()
252 original_photo_calib = input_exposure.getPhotoCalib()
253 original_wcs = input_exposure.getWcs()
254 input_exposure.setPsf(psf)
255 input_exposure.setPhotoCalib(photo_calib)
256 input_exposure.setWcs(wcs)
257 bbox = input_exposure.getBBox()
258 photo_calib = input_exposure.getPhotoCalib()
259 elif isinstance(input_exposure, CellCoadd):
260 if wcs is None:
261 wcs = input_exposure.sky_projection.to_legacy()
262 bbox = input_exposure.bbox.to_legacy()
263 photo_calib = BaseField.make_legacy_photo_calib(input_exposure.unit)
264 else:
265 raise ValueError(f"Unsupported {type(input_exposure)=} for injection")
267 # Make empty table if none supplied to support process_all_data_ids.
268 if len(injection_catalogs) == 0:
269 if self.config.process_all_data_ids:
270 injection_catalogs = [Table(names=["ra", "dec", "source_type"])]
271 else:
272 corners = [Point2D(x) for x in input_exposure.getBBox().getCorners()]
273 coords = wcs.pixelToSky(corners)
274 ras = [x.getRa().asDegrees() for x in coords]
275 decs = [x.getDec().asDegrees() for x in coords]
276 raise RuntimeError(
277 "No injection sources overlap the region bounded by "
278 f"{np.min(ras):.2f} <= RA < {np.max(ras):.2f}, "
279 f"{np.min(decs):.2f} <= Dec < {np.max(decs):.2f} (degrees). "
280 "Check injection catalog coverage."
281 )
283 # Consolidate injection catalogs and compose main injection catalog.
284 injection_catalog = self._compose_injection_catalog(injection_catalogs)
286 # Mapping between standard column names and configured names/units.
287 column_mapping = {
288 "ra": (self.config.col_ra, units.deg),
289 "dec": (self.config.col_dec, units.deg),
290 "source_type": (self.config.col_source_type, None),
291 "mag": (self.config.col_mag, units.mag),
292 "stamp": (self.config.col_stamp, None),
293 "draw_size": (self.config.col_draw_size, units.pix),
294 "trail_length": (self.config.col_trail_length, units.pix),
295 }
297 # Standardize injection catalog column names and units.
298 injection_catalog = self._standardize_columns(
299 injection_catalog,
300 column_mapping,
301 wcs.getPixelScale(bbox.getCenter()).asArcseconds(),
302 )
304 # Clean the injection catalog of sources which are not injectable.
305 injection_catalog = self._clean_sources(injection_catalog, input_exposure, bbox=bbox, wcs=wcs)
307 # Injection binary flag lookup dictionary.
308 binary_flags = {
309 "MAG_BAD": 0,
310 "TYPE_UNKNOWN": 1,
311 "SERSIC_EXTREME": 2,
312 "NO_OVERLAP": 3,
313 "FFT_SIZE_ERROR": 4,
314 "PSF_COMPUTE_ERROR": 5,
315 }
317 # Check that sources in the injection catalog are able to be injected.
318 injection_catalog = self._check_sources(injection_catalog, binary_flags)
320 # Inject sources into input_exposure.
321 good_injections: list[bool] = injection_catalog["injection_flag"] == 0
322 good_injections_index = [i for i, val in enumerate(good_injections) if val]
323 num_injection_sources = np.sum(good_injections)
324 if num_injection_sources > 0:
325 object_generator = generate_galsim_objects(
326 injection_catalog=injection_catalog[good_injections],
327 photo_calib=photo_calib,
328 wcs=wcs,
329 fits_alignment=self.config.fits_alignment,
330 stamp_prefix=self.config.stamp_prefix,
331 logger=self.log,
332 )
333 (
334 draw_sizes,
335 common_bounds,
336 fft_size_errors,
337 psf_compute_errors,
338 ) = inject_galsim_objects_into_exposure(
339 input_exposure,
340 object_generator,
341 mask_plane_name=self.config.mask_plane_name,
342 calib_flux_radius=self.config.calib_flux_radius,
343 draw_size_max=10000, # TODO: replace draw_size logic with GS logic.
344 inject_variance=self.config.inject_variance,
345 add_noise=self.config.add_noise,
346 noise_seed=self.config.noise_seed,
347 bad_mask_names=list(self.config.bad_mask_names),
348 logger=self.log,
349 )
350 # Add inject_galsim_objects_into_exposure outputs into output cat.
351 common_areas = [x.area() if x is not None else None for x in common_bounds]
352 for i, (draw_size, common_area, fft_size_error, psf_compute_error) in enumerate(
353 zip(draw_sizes, common_areas, fft_size_errors, psf_compute_errors)
354 ):
355 injection_catalog["injection_draw_size"][good_injections_index[i]] = draw_size
356 if common_area == 0:
357 injection_catalog["injection_flag"][good_injections_index[i]] += (
358 2 ** binary_flags["NO_OVERLAP"]
359 )
360 if fft_size_error:
361 injection_catalog["injection_flag"][good_injections_index[i]] += (
362 2 ** binary_flags["FFT_SIZE_ERROR"]
363 )
364 if psf_compute_error:
365 injection_catalog["injection_flag"][good_injections_index[i]] += (
366 2 ** binary_flags["PSF_COMPUTE_ERROR"]
367 )
368 num_injected_sources = np.sum(injection_catalog["injection_flag"] == 0)
369 num_skipped_sources = np.sum(injection_catalog["injection_flag"] != 0)
370 grammar1 = "source" if num_injection_sources == 1 else "sources"
371 grammar2 = "source" if num_skipped_sources == 1 else "sources"
373 injection_flags = np.array(injection_catalog["injection_flag"])
374 num_injection_flags = [np.sum((injection_flags & 2**x) > 0) for x in binary_flags.values()]
375 if np.sum(num_injection_flags) > 0:
376 injection_flag_report = ": " + ", ".join(
377 [f"{x}({y})" for x, y in zip(binary_flags.keys(), num_injection_flags) if y > 0]
378 )
379 else:
380 injection_flag_report = ""
381 self.log.info(
382 "Injected %d of %d potential %s. %d %s flagged and skipped%s.",
383 num_injected_sources,
384 num_injection_sources,
385 grammar1,
386 num_skipped_sources,
387 grammar2,
388 injection_flag_report,
389 )
390 elif num_injection_sources == 0 and self.config.process_all_data_ids:
391 self.log.warning("No sources to be injected for this DatasetRef; processing anyway.")
392 input_exposure.mask.addMaskPlane(self.config.mask_plane_name)
393 mask_plane_core_name = self.config.mask_plane_name + "_CORE"
394 input_exposure.mask.addMaskPlane(mask_plane_core_name)
395 self.log.info(
396 "Adding %s and %s mask planes to the exposure.",
397 self.config.mask_plane_name,
398 mask_plane_core_name,
399 )
400 else:
401 raise RuntimeError(
402 "No sources to be injected for this DatasetRef, and process_all_data_ids is False."
403 )
405 # Restore original input_exposure calibrated data.
406 if is_exposure:
407 input_exposure.setPsf(original_psf)
408 input_exposure.setPhotoCalib(original_photo_calib)
409 input_exposure.setWcs(original_wcs)
411 # Add injection provenance and injection flags metadata.
412 metadata = input_exposure.getMetadata()
413 input_dataset_type = self.config.connections.input_exposure.format(
414 **self.config.connections.toDict()
415 )
416 metadata.set("INJECTED", input_dataset_type, "Initial source injection dataset type")
417 input_exposure.getInfo().setVisitInfo(input_exposure.visitInfo.copyWith(hasSimulatedContent=True))
418 for flag, value in sorted(binary_flags.items(), key=lambda item: item[1]):
419 injection_catalog.meta[flag] = value
421 output_struct = Struct(output_exposure=input_exposure, output_catalog=injection_catalog)
422 return output_struct
424 def _compose_injection_catalog(self, injection_catalogs):
425 """Consolidate injection catalogs and compose main injection catalog.
427 If multiple injection catalogs are input, all catalogs are
428 concatenated together.
430 A running injection_id, specific to this dataset ref, is assigned to
431 each source in the output injection catalog if not provided.
433 Parameters
434 ----------
435 injection_catalogs : `list` [`astropy.table.Table`]
436 Set of synthetic source catalogs to concatenate.
438 Returns
439 -------
440 injection_catalog : `astropy.table.Table`
441 Catalog of sources to be injected.
442 """
443 self.config = cast(BaseInjectConfig, self.config)
445 # Generate injection IDs (if not provided) and injection flag column.
446 injection_data = vstack(injection_catalogs, metadata_conflicts="silent")
447 if "injection_id" in injection_data.columns:
448 injection_id = injection_data["injection_id"]
449 injection_data.remove_column("injection_id")
450 else:
451 injection_id = range(len(injection_data))
452 injection_header = Table(
453 {
454 "injection_id": injection_id,
455 "injection_flag": np.zeros(len(injection_data), dtype=int),
456 "injection_draw_size": np.zeros(len(injection_data), dtype=int),
457 }
458 )
460 # Construct final injection catalog.
461 injection_catalog = hstack([injection_header, injection_data])
462 injection_catalog["source_type"] = injection_catalog["source_type"].astype(str)
464 # Log and return.
465 num_injection_catalogs = np.sum([len(table) > 0 for table in injection_catalogs])
466 grammar1 = "source" if len(injection_catalog) == 1 else "sources"
467 grammar2 = "trixel" if num_injection_catalogs == 1 else "trixels"
468 self.log.info(
469 "Retrieved %d injection %s from %d HTM %s.",
470 len(injection_catalog),
471 grammar1,
472 num_injection_catalogs,
473 grammar2,
474 )
475 return injection_catalog
477 def _standardize_columns(self, injection_catalog, column_mapping, pixel_scale):
478 """Standardize injection catalog column names and units.
480 Use config variables to standardize the expected columns and column
481 names in the input injection catalog. This method replaces all core
482 column names in the config with hard-coded internal names.
484 Only a core group of column names are standardized; additional column
485 names will not be modified. If certain parameters are needed (i.e.,
486 by GalSim), these columns must be given exactly as required in the
487 appropriate units. Refer to the configuration documentation for more
488 details.
490 Parameters
491 ----------
492 injection_catalog : `astropy.table.Table`
493 A catalog of sources to be injected.
494 column_mapping : `dict` [`str`, `tuple` [`str`, `astropy.units.Unit`]]
495 A dictionary mapping standard column names to the configured column
496 names and units.
497 pixel_scale : `float`
498 Pixel scale of the exposure in arcseconds per pixel.
500 Returns
501 -------
502 injection_catalog : `astropy.table.Table`
503 The standardized catalog of sources to be injected.
504 """
505 self.config = cast(BaseInjectConfig, self.config)
507 pixel_scale_equivalency = units.pixel_scale(
508 Quantity(pixel_scale, units.arcsec / units.pix) # type: ignore
509 )
510 for standard_col, (configured_col, unit) in column_mapping.items():
511 # Rename columns if necessary.
512 if configured_col in injection_catalog.colnames:
513 injection_catalog.rename_column(configured_col, standard_col)
514 # Attempt to convert to our desired units, then remove units.
515 if standard_col in injection_catalog.columns and unit:
516 try:
517 injection_catalog[standard_col] = (
518 injection_catalog[standard_col].to(unit, pixel_scale_equivalency).value
519 )
520 except UnitConversionError:
521 pass
522 return Table(injection_catalog)
524 def _clean_sources(self, injection_catalog, input_exposure, bbox=None, wcs=None):
525 """Clean the injection catalog of sources which are not injectable.
527 This method will remove sources which are not injectable for a variety
528 of reasons, namely: sources which fall outside the padded exposure
529 bounding box or sources not selected by virtue of their evaluated
530 selection criteria.
532 If the input injection catalog contains x/y inputs but does not contain
533 RA/Dec inputs, WCS information will be used to generate RA/Dec sky
534 coordinate information and appended to the injection catalog.
536 Parameters
537 ----------
538 injection_catalog : `astropy.table.Table`
539 The catalog of sources to be injected.
540 input_exposure : `lsst.afw.image.ExposureF`
541 The exposure to inject sources into.
543 Returns
544 -------
545 injection_catalog : `astropy.table.Table`
546 Updated injection catalog containing *x* and *y* pixel coordinates,
547 and cleaned to only include injection sources which fall within the
548 bounding box of the input exposure dilated by *trim_padding*.
549 """
550 self.config = cast(BaseInjectConfig, self.config)
552 if wcs is None:
553 wcs = input_exposure.getWcs()
555 if bbox is None:
556 bbox = input_exposure.getBBox()
558 # Exit early if there are no sources to inject.
559 if len(injection_catalog) == 0:
560 self.log.info("Catalog cleaning not applied to empty injection catalog.")
561 return injection_catalog
563 sources_to_keep = np.ones(len(injection_catalog), dtype=bool)
565 # Determine centroids and remove sources outside the padded bbox.
566 has_sky = {"ra", "dec"} <= set(injection_catalog.columns)
567 has_pixel = {"x", "y"} <= set(injection_catalog.columns)
568 # Input catalog must contain either RA/Dec OR x/y.
569 # If only x/y given, RA/Dec will be calculated.
570 if not has_sky and has_pixel:
571 begin_x, begin_y = bbox.getBegin()
572 ras, decs = wcs.pixelToSkyArray(
573 begin_x + injection_catalog["x"].astype(float),
574 begin_y + injection_catalog["y"].astype(float),
575 degrees=True,
576 )
577 injection_catalog["ra"] = ras
578 injection_catalog["dec"] = decs
579 injection_catalog["x"] += begin_x
580 injection_catalog["y"] += begin_y
581 has_sky = True
582 elif not has_sky and not has_pixel:
583 self.log.warning("No spatial coordinates found in injection catalog; cannot inject any sources!")
584 if has_sky:
585 if self.config.trim_padding:
586 bbox.grow(int(self.config.trim_padding))
587 is_contained = bbox_contains_sky_coords(
588 bbox, wcs, injection_catalog["ra"] * units.deg, injection_catalog["dec"] * units.deg
589 )
590 sources_to_keep &= is_contained
591 if (num_not_contained := np.sum(~is_contained)) > 0:
592 grammar = ("source", "a centroid") if num_not_contained == 1 else ("sources", "centroids")
593 self.log.info(
594 "Identified %d injection %s with %s outside the padded image bounding box.",
595 num_not_contained,
596 grammar[0],
597 grammar[1],
598 )
600 # Remove sources by boolean selection flag.
601 if self.config.selection:
602 visit = input_exposure.getInfo().getVisitInfo().getId()
603 selected = eval(self.config.selection.format(visit=visit))
604 sources_to_keep &= selected
605 if (num_not_selected := np.sum(~selected)) >= 0:
606 grammar = ["source", "was"] if num_not_selected == 1 else ["sources", "were"]
607 self.log.warning(
608 "Identified %d injection %s that %s not selected.",
609 num_not_selected,
610 grammar[0],
611 grammar[1],
612 )
614 # Print final cleaning report and return.
615 num_cleaned_total = np.sum(~sources_to_keep)
616 grammar = "source" if len(sources_to_keep) == 1 else "sources"
617 self.log.info(
618 "Catalog cleaning removed %d of %d %s; %d remaining for catalog checking.",
619 num_cleaned_total,
620 len(sources_to_keep),
621 grammar,
622 np.sum(sources_to_keep),
623 )
624 injection_catalog = injection_catalog[sources_to_keep]
625 return injection_catalog
627 def _check_sources(self, injection_catalog, binary_flags):
628 """Check that sources in the injection catalog are able to be injected.
630 This method will check that sources in the injection catalog are able
631 to be injected, and will flag them if not. Checks will be made on a
632 number of parameters, including magnitude, source type and Sérsic index
633 (where relevant).
635 Legacy profile types will be renamed to their standardized GalSim
636 equivalents; any source profile types that are not GalSim classes will
637 be flagged.
639 Note: Unlike the cleaning method, no sources are actually removed here.
640 Instead, a binary flag is set in the *injection_flag* column for each
641 source. Only unflagged sources will be generated for source injection.
643 Parameters
644 ----------
645 injection_catalog : `astropy.table.Table`
646 Catalog of sources to be injected.
647 binary_flags : `dict` [`str`, `int`]
648 Dictionary of binary flags to be used in the injection_flag column.
650 Returns
651 -------
652 injection_catalog : `astropy.table.Table`
653 The cleaned catalog of sources to be injected.
654 """
655 self.config = cast(BaseInjectConfig, self.config)
657 # Exit early if there are no sources to inject.
658 if len(injection_catalog) == 0:
659 self.log.info("Catalog checking not applied to empty injection catalog.")
660 return injection_catalog
662 # Flag erroneous magnitude values (missing mag data or NaN mag values).
663 if "mag" not in injection_catalog.columns:
664 # Check injection_catalog has a mag column.
665 self.log.warning("No magnitude data found in injection catalog; cannot inject any sources!")
666 injection_catalog["injection_flag"] += 2 ** binary_flags["MAG_BAD"]
667 else:
668 # Check that all input mag values are finite.
669 mag_array = np.isfinite(ma.array(injection_catalog["mag"]))
670 bad_mag = ~(mag_array.data * ~mag_array.mask)
671 if (num_bad_mag := np.sum(bad_mag)) > 0:
672 grammar = "source" if num_bad_mag == 1 else "sources"
673 self.log.warning(
674 "Flagging %d injection %s that do not have a finite magnitude.", num_bad_mag, grammar
675 )
676 injection_catalog["injection_flag"][bad_mag] += 2 ** binary_flags["MAG_BAD"]
678 # Replace legacy source types with standardized profile names.
679 injection_catalog["source_type"] = injection_catalog["source_type"].astype("O")
680 replace_dict = {"Star": "DeltaFunction"}
681 for legacy_type, standard_type in replace_dict.items():
682 legacy_matches = injection_catalog["source_type"] == legacy_type
683 if np.any(legacy_matches):
684 injection_catalog["source_type"][legacy_matches] = standard_type
685 injection_catalog["source_type"] = injection_catalog["source_type"].astype(str)
687 # Flag source types not supported by GalSim.
688 input_source_types = set(injection_catalog["source_type"])
689 for input_source_type in input_source_types:
690 if input_source_type not in _ALLOWED_SOURCE_TYPES:
691 unknown_source_types = injection_catalog["source_type"] == input_source_type
692 grammar = "source" if np.sum(unknown_source_types) == 1 else "sources"
693 self.log.warning(
694 "Flagging %d injection %s with an unsupported source type: %s.",
695 np.sum(unknown_source_types),
696 grammar,
697 input_source_type,
698 )
699 injection_catalog["injection_flag"][unknown_source_types] += 2 ** binary_flags["TYPE_UNKNOWN"]
701 # Flag extreme Sersic index sources.
702 if "n" in injection_catalog.columns:
703 min_n = galsim.Sersic._minimum_n
704 max_n = galsim.Sersic._maximum_n
705 n_vals = injection_catalog["n"]
706 extreme_sersics = (n_vals <= min_n) | (n_vals >= max_n)
707 if (num_extreme_sersics := np.sum(extreme_sersics)) > 0:
708 grammar = "source" if num_extreme_sersics == 1 else "sources"
709 self.log.warning(
710 "Flagging %d injection %s with a Sersic index outside the range %.1f <= n <= %.1f.",
711 num_extreme_sersics,
712 grammar,
713 min_n,
714 max_n,
715 )
716 injection_catalog["injection_flag"][extreme_sersics] += 2 ** binary_flags["SERSIC_EXTREME"]
718 # Print final cleaning report.
719 num_flagged_total = np.sum(injection_catalog["injection_flag"] != 0)
720 grammar = "source" if len(injection_catalog) == 1 else "sources"
721 self.log.info(
722 "Catalog checking flagged %d of %d %s; %d remaining for source generation.",
723 num_flagged_total,
724 len(injection_catalog),
725 grammar,
726 np.sum(injection_catalog["injection_flag"] == 0),
727 )
728 return injection_catalog