Coverage for python/lsst/source/injection/inject

1# This file is part of source_injection.

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.

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__ = ["generate_galsim_objects", "inject_galsim_objects_into_exposure"]

26import os

27from collections import Counter

28from collections.abc import Generator

29from typing import Any

31import galsim

32import numpy as np

33import numpy.ma as ma

34from astropy.io import fits

35from astropy.table import Table

36from galsim import GalSimFFTSizeError

38from lsst.afw.geom import SkyWcs

39from lsst.afw.image import ExposureF, PhotoCalib

40from lsst.geom import Box2I, Point2D, Point2I, SpherePoint, arcseconds, degrees

41from lsst.pex.exceptions import InvalidParameterError, LogicError

44def get_object_data(source_data: dict[str, Any], object_class: galsim.GSObject) -> dict[str, Any]:

45 """Assemble a dictionary of allowed keyword arguments and their

46 corresponding values to use when constructing a GSObject.

48 Parameters

49 ----------

50 source_data : `dict` [`str`, `Any`]

51 Dictionary of source data.

52 object_class : `galsim.gsobject.GSObject`

53 Class of GSObject to match against.

55 Returns

56 -------

57 object_data : `dict` [`str`, `Any`]

58 Dictionary of source data to pass to the GSObject constructor.

59 """

60 req = getattr(object_class, "_req_params", {})

61 opt = getattr(object_class, "_opt_params", {})

62 single = getattr(object_class, "_single_params", {})

64 object_data = {}

66 # Check required args.

67 for key in req:

68 if key not in source_data:

69 raise ValueError(f"Required parameter {key} not found in input catalog.")

70 object_data[key] = source_data[key]

72 # Optional args.

73 for key in opt:

74 if key in source_data:

75 object_data[key] = source_data[key]

77 # Single args.

78 for s in single:

79 count = 0

80 for key in s:

81 if key in source_data:

82 count += 1

83 if count > 1:

84 raise ValueError(f"Only one of {s.keys()} allowed for type {object_class}.")

85 object_data[key] = source_data[key]

86 if count == 0:

87 raise ValueError(f"One of the args {s.keys()} is required for type {object_class}.")

89 return object_data

92def get_shear_data(

93 source_data: dict[str, Any],

94 shear_attributes: list[str] = [

95 "g1",

96 "g2",

97 "g",

98 "e1",

99 "e2",

100 "e",

101 "eta1",

102 "eta2",

103 "eta",

104 "q",

105 "beta",

106 "shear",

107 ],

108) -> dict[str, Any]:

109 """Assemble a dictionary of allowed keyword arguments and their

110 corresponding values to use when constructing a Shear.

111

112 Parameters

113 ----------

114 source_data : `dict` [`str`, `Any`]

115 Dictionary of source data.

116 shear_attributes : `list` [`str`], optional

117 List of allowed shear attributes.

118

119 Returns

120 -------

121 shear_data : `dict` [`str`, `Any`]

122 Dictionary of source data to pass to the Shear constructor.

123 """

124 shear_params = set(shear_attributes) & set(source_data.keys())

125 shear_data = {}

126 for shear_param in shear_params:

127 if shear_param == "beta":

128 shear_data.update({shear_param: source_data[shear_param] * galsim.degrees})

129 else:

130 shear_data.update({shear_param: source_data[shear_param]})

131 return shear_data

132

133

134def generate_galsim_objects(

135 injection_catalog: Table,

136 photo_calib: PhotoCalib,

137 wcs: SkyWcs,

138 fits_alignment: str,

139 stamp_prefix: str = "",

140 logger: Any | None = None,

141) -> Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None]:

142 """Generate GalSim objects from an injection catalog.

143

144 Parameters

145 ----------

146 injection_catalog : `astropy.table.Table`

147 Table of sources to be injected.

148 photo_calib : `lsst.afw.image.PhotoCalib`

149 Photometric calibration used to calibrate injected sources.

150 wcs : `lsst.afw.geom.SkyWcs`

151 WCS used to calibrate injected sources.

152 fits_alignment : `str`

153 Alignment of the FITS image to the WCS. Allowed values: "wcs", "pixel".

154 stamp_prefix : `str`

155 Prefix to add to the stamp name.

156 logger : `lsst.utils.logging.LsstLogAdapter`, optional

157 Logger to use for logging messages.

158

159 Yields

160 ------

161 sky_coords : `lsst.geom.SpherePoint`

162 RA/Dec coordinates of the source.

163 pixel_coords : `lsst.geom.Point2D`

164 Pixel coordinates of the source.

165 draw_size : `int`

166 Size of the stamp to draw.

167 object : `galsim.gsobject.GSObject`

168 A fully specified and transformed GalSim object.

169 """

170 if logger:

171 source_types = Counter(injection_catalog["source_type"]) # type: ignore

172 grammar0 = "source" if len(injection_catalog) == 1 else "sources"

173 grammar1 = "type" if len(source_types) == 1 else "types"

174 logger.info(

175 "Generating %d injection %s consisting of %d unique %s: %s.",

176 len(injection_catalog),

177 grammar0,

178 len(source_types),

179 grammar1,

180 ", ".join(f"{k}({v})" for k, v in source_types.items()),

181 )

182 for source_data_full in injection_catalog:

183 items = dict(source_data_full).items() # type: ignore

184 source_data = {k: v for (k, v) in items if v is not ma.masked}

185 try:

186 sky_coords = SpherePoint(float(source_data["ra"]), float(source_data["dec"]), degrees)

187 except KeyError:

188 sky_coords = wcs.pixelToSky(float(source_data["x"]), float(source_data["y"]))

189 try:

190 pixel_coords = Point2D(source_data["x"], source_data["y"])

191 except KeyError:

192 pixel_coords = wcs.skyToPixel(sky_coords)

193 try:

194 inst_flux = photo_calib.magnitudeToInstFlux(source_data["mag"], pixel_coords)

195 except LogicError:

196 continue

197 try:

198 draw_size = int(source_data["draw_size"])

199 except KeyError:

200 draw_size = 0

201

202 if source_data["source_type"] == "Stamp":

203 stamp_file = stamp_prefix + source_data["stamp"]

204 object = make_galsim_stamp(stamp_file, fits_alignment, wcs, sky_coords, inst_flux)

205 elif source_data["source_type"] == "Trail":

206 object = make_galsim_trail(source_data, wcs, sky_coords, inst_flux)

207 else:

208 object = make_galsim_object(source_data, source_data["source_type"], inst_flux, logger)

209

210 yield sky_coords, pixel_coords, draw_size, object

211

212

213def make_galsim_object(

214 source_data: dict[str, Any],

215 source_type: str,

216 inst_flux: float,

217 logger: Any | None = None,

218) -> galsim.gsobject.GSObject:

219 """Make a generic GalSim object from a collection of source data.

220

221 Parameters

222 ----------

223 source_data : `dict` [`str`, `Any`]

224 Dictionary of source data.

225 source_type : `str`

226 Type of the source, corresponding to a GalSim class.

227 inst_flux : `float`

228 Instrumental flux of the source.

229 logger : `~lsst.utils.logging.LsstLogAdapter`, optional

230

231 Returns

232 -------

233 object : `galsim.gsobject.GSObject`

234 A fully specified and transformed GalSim object.

235 """

236 # Populate the non-shear and non-flux parameters.

237 object_class = getattr(galsim, source_type)

238 object_data = get_object_data(source_data, object_class)

239 object = object_class(**object_data)

240 # Create a version of the object with an area-preserving shear applied.

241 shear_data = get_shear_data(source_data)

242 if shear_data:

243 try:

244 object = object.shear(**shear_data)

245 except TypeError as err:

246 if logger:

247 logger.warning("Cannot apply shear to GalSim object: %s", err)

248 pass

249 # Apply the instrumental flux and return.

250 object = object.withFlux(inst_flux)

251 return object

252

253

254def make_galsim_trail(

255 source_data: dict[str, Any],

256 wcs: SkyWcs,

257 sky_coords: SpherePoint,

258 inst_flux: float,

259 trail_thickness: float = 1e-6,

260) -> galsim.gsobject.GSObject:

261 """Make a trail with GalSim from a collection of source data.

262

263 Parameters

264 ----------

265 source_data : `dict` [`str`, `Any`]

266 Dictionary of source data.

267 wcs : `lsst.afw.geom.SkyWcs`

268 World coordinate system.

269 sky_coords : `lsst.geom.SpherePoint`

270 Sky coordinates of the source.

271 inst_flux : `float`

272 Instrumental flux of the source.

273 trail_thickness : `float`

274 Thickness of the trail in pixels.

275

276 Returns

277 -------

278 object : `galsim.gsobject.GSObject`

279 A fully specified and transformed GalSim object.

280 """

281 # Make a 'thin' box to mimic a line surface brightness profile of default

282 # thickness = 1e-6 (i.e., much thinner than a pixel)

283 object = galsim.Box(source_data["trail_length"], trail_thickness)

284 try:

285 object = object.rotate(source_data["beta"] * galsim.degrees)

286 except KeyError:

287 pass

288 object = object.withFlux(inst_flux * source_data["trail_length"]) # type: ignore

289 # GalSim objects are assumed to be in sky-coords. As we want the trail to

290 # appear as defined above in image-coords, we must transform it here.

291 linear_wcs = wcs.linearizePixelToSky(sky_coords, arcseconds)

292 mat = linear_wcs.getMatrix()

293 object = object.transform(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1])

294 # Rescale flux; transformation preserves surface brightness, not flux.

295 object *= 1.0 / np.abs(mat[0, 0] * mat[1, 1] - mat[0, 1] * mat[1, 0])

296 return object # type: ignore

297

298

299def make_galsim_stamp(

300 stamp_file: str,

301 fits_alignment: str,

302 wcs: SkyWcs,

303 sky_coords: SpherePoint,

304 inst_flux: float,

305) -> galsim.gsobject.GSObject:

306 """Make a postage stamp with GalSim from a FITS file.

307

308 Parameters

309 ----------

310 stamp_file : `str`

311 Path to the FITS file containing the postage stamp.

312 fits_alignment : `str`

313 Alignment of the FITS image to the WCS. Allowed values: "wcs", "pixel".

314 wcs : `lsst.afw.geom.SkyWcs`

315 World coordinate system.

316 sky_coords : `lsst.geom.SpherePoint`

317 Sky coordinates of the source.

318 inst_flux : `float`

319 Instrumental flux of the source.

320

321 Returns

322 -------

323 object: `galsim.gsobject.GSObject`

324 A fully specified and transformed GalSim object.

325 """

326 stamp_file = stamp_file.strip()

327 if os.path.exists(stamp_file):

328 with fits.open(stamp_file) as hdul:

329 hdu_images = [hdu.is_image and hdu.size > 0 for hdu in hdul] # type: ignore

330 if any(hdu_images):

331 stamp_data = galsim.fits.read(stamp_file, read_header=True, hdu=np.where(hdu_images)[0][0])

332 else:

333 raise RuntimeError(f"Cannot find image in input FITS file {stamp_file}.")

334 match fits_alignment:

335 case "wcs":

336 # galsim.fits.read will always attach a WCS to its output.

337 # If it can't find a WCS in the FITS header, then it

338 # defaults to scale = 1.0 arcsec / pix. If that's the scale

339 # then we need to check if it was explicitly set or if it's

340 # just the default. If it's just the default then we should

341 # raise an exception.

342 if is_wcs_galsim_default(stamp_data.wcs, stamp_data.header): # type: ignore

343 raise RuntimeError(f"Cannot find WCS in input FITS file {stamp_file}")

344 case "pixel":

345 # We need to set stamp_data.wcs to the local target WCS.

346 linear_wcs = wcs.linearizePixelToSky(sky_coords, arcseconds)

347 mat = linear_wcs.getMatrix()

348 stamp_data.wcs = galsim.JacobianWCS( # type: ignore

349 mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]

350 )

351 object = galsim.InterpolatedImage(stamp_data, calculate_stepk=False)

352 object = object.withFlux(inst_flux)

353 return object # type: ignore

354 else:

355 raise RuntimeError(f"Cannot locate input FITS postage stamp {stamp_file}.")

356

357

358def is_wcs_galsim_default(

359 wcs: galsim.fitswcs.GSFitsWCS,

360 hdr: galsim.fits.FitsHeader,

361) -> bool:

362 """Decide if wcs = galsim.PixelScale(1.0) is explicitly present in header,

363 or if it's just the GalSim default.

364

365 Parameters

366 ----------

367 wcs : galsim.fitswcs.GSFitsWCS

368 Potentially default WCS.

369 hdr : galsim.fits.FitsHeader

370 Header as read in by GalSim.

371

372 Returns

373 -------

374 is_default : bool

375 True if default, False if explicitly set in header.

376 """

377 if wcs != galsim.PixelScale(1.0):

378 return False

379 if hdr.get("GS_WCS") is not None:

380 return False

381 if hdr.get("CTYPE1", "LINEAR") == "LINEAR":

382 return not any(k in hdr for k in ["CD1_1", "CDELT1"])

383 for wcs_type in galsim.fitswcs.fits_wcs_types:

384 # If one of these succeeds, then assume result is explicit.

385 try:

386 wcs_type._readHeader(hdr)

387 return False

388 except Exception:

389 pass

390 else:

391 return not any(k in hdr for k in ["CD1_1", "CDELT1"])

392

393

394def get_gains_from_amps(

395 exposure: ExposureF,

396) -> tuple[list[Box2I], list[float]] | tuple[None, None]:

397 """Retrieve the bounding boxes and gains for each amplifier, if applicable.

398 If this cannot be done, return None.

399

400 Parameters

401 ----------

402 exposure : `lsst.afw.image.ExposureF`

403 The exposure to inject synthetic sources into.

404

405 Returns

406 -------

407 bboxes : `list` [`lsst.geom.Box2I`]

408 Bounding boxes for each amplifier region in the exposure.

409 gains : `list` [`float`]

410 The gain recorded for each amplifier in this exposure.

411 """

412 try:

413 amplifiers = exposure.getDetector().getAmplifiers()

414 bboxes = [amplifier.getBBox() for amplifier in amplifiers]

415 gains = [amplifier.getGain() for amplifier in amplifiers]

416 except AttributeError:

417 return None, None

418 return bboxes, gains

419

420

421def infer_gain_from_image(

422 exposure: ExposureF,

423 bad_mask_names: list[str] | None = None,

424) -> float:

425 """Fit a straight line to image flux vs. estimated variance in each pixel

426 of an exposure, and from that infer an estimated average gain.

427

428 Parameters

429 ----------

430 exposure : `lsst.afw.image.ExposureF`

431 The exposure to inject synthetic sources into.

432 bad_mask_names : `list` [`str`], optional

433 Names of any mask plane bits that could indicate a pathological flux

434 vs. variance relationship.

435

436 Returns

437 -------

438 gain : `float`

439 The estimated average gain across all pixels.

440 """

441 image_arr = exposure.image.array

442 var_arr = exposure.variance.array

443 mask_arr = exposure.mask.array

444 good = np.isfinite(image_arr) & np.isfinite(var_arr)

445 if bad_mask_names is not None:

446 bad_mask_bit_mask = exposure.mask.getPlaneBitMask(bad_mask_names)

447 good &= (mask_arr.astype(int) & bad_mask_bit_mask) == 0

448 fit = np.polyfit(image_arr[good], var_arr[good], deg=1)

449 return 1.0 / fit[0]

450

451

452def get_gain_map(

453 exposure: ExposureF,

454 full_bounds: galsim.BoundsI,

455 bad_mask_names: list[str] | None = None,

456) -> tuple[galsim.Image, bool]:

457 """Retrieve gains for each pixel.

458

459 Parameters

460 ----------

461 exposure : `lsst.afw.image.ExposureF`

462 The exposure to inject synthetic sources into.

463 full_bounds : `galsim.BoundsI`

464 Bounding box covering the entire exposure.

465 bad_mask_names : `list` [`str`], optional

466 Names of any mask plane bits that could indicate a pathological flux

467 vs. variance relationship.

468

469 Returns

470 -------

471 gain_map : `galsim.Image`

472 The gain to use in each pixel.

473 is_single_CCD : `bool`

474 Whether the exposure represents a single CCD, having well-defined

475 amplifier regions with recorded gains.

476 """

477 gain_map = galsim.Image(full_bounds, dtype=float)

478 is_single_CCD = True

479

480 bboxes, gains = get_gains_from_amps(exposure)

481 if bboxes is None:

482 # Separate amplifier bounding boxes not found.

483 # This exposure does not represent a single CCD.

484 is_single_CCD = False

485 bboxes = [full_bounds]

486 gains = [infer_gain_from_image(exposure, bad_mask_names)]

487

488 assert gains is not None # Make mypy happy.

489 for bbox, gain in zip(bboxes, gains):

490 if is_single_CCD:

491 bounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY)

492 else:

493 bounds = galsim.BoundsI(bbox.getXMin(), bbox.getXMax(), bbox.getYMin(), bbox.getYMax())

494 gain_map[bounds].array[:] = 1.0 / gain

495 return gain_map, is_single_CCD

496

497

498def add_noise_to_galsim_image(

499 image_template: galsim.Image,

500 var_template: galsim.Image,

501 is_single_CCD: bool,

502 gain_map: galsim.Image,

503 object_common_bounds: galsim.BoundsI,

504 noise_seed: int | None = None,

505 logger: Any | None = None,

506) -> tuple[galsim.Image, galsim.Image]:

507 """Add noise to a supplied image, and adjust the estimated variance

508 accordingly. If the image represents a single CCD exposure, add Poisson

509 shot noise. Otherwise add Gaussian noise.

510

511 Parameters

512 ----------

513 image_template : `galsim.Image`

514 The image to add noise to.

515 var_template : `galsim.Image`

516 The variance to use for the noise generating process in each pixel.

517 is_single_CCD : `bool`

518 Whether the exposure represents a single CCD, having well-defined

519 amplifier regions with recorded gains.

520 gain_map : `galsim.Image`

521 The gain to use in each pixel of the full exposure.

522 object_common_bounds : `galsim.BoundsI`

523 Specific region of the full exposure to which this image belongs.

524 noise_seed : `int`, optional

525 Random number generator seed to use for the noise realization.

526 logger : `lsst.utils.logging.LsstLogAdapter`, optional

527 Logger to use for logging messages.

528

529 Returns

530 -------

531 image_template : `galsim.Image`

532 The image with noise added.

533 var_template : `galsim.Image`

534 The estimated variance, adjusted in response to the added noise.

535 """

536 noise_template = var_template.copy()

537 negative_variance = var_template.array < 0

538 nonfinite_variance = ~np.isfinite(var_template.array)

539 if np.any(negative_variance):

540 if logger:

541 logger.debug("Setting negative-variance pixels to 0 variance for noise generation.")

542 noise_template.array[negative_variance] = 0

543 if np.any(nonfinite_variance):

544 if logger:

545 logger.debug("Setting non-finite-variance pixels to 0 variance for noise generation.")

546 noise_template.array[nonfinite_variance] = 0

547

548 rng = galsim.BaseDeviate(noise_seed)

549

550 if is_single_CCD:

551 # Treat noise_template, which is the expected amount of

552 # injected flux divided by the gain, as the Poisson mean of

553 # the number of photons collected by each pixel.

554 noise = galsim.PoissonNoise(rng)

555 noise_template.addNoise(noise)

556 # Scale the photon count by the gain to get the amount of

557 # image flux to inject.

558 image_template = noise_template.copy()

559 image_template /= gain_map[object_common_bounds]

560 # Restore the original variance values

561 # of the bad-variance pixels.

562 bad_var = negative_variance | nonfinite_variance

563 noise_template.array[bad_var] = var_template.array[bad_var]

564 var_template = noise_template

565 else:

566 noise = galsim.VariableGaussianNoise(rng, noise_template)

567 image_template.addNoise(noise)

568 var_template = image_template.copy()

569 var_template *= gain_map[object_common_bounds]

570

571 return image_template, var_template

572

573

574def inject_galsim_objects_into_exposure(

575 exposure: ExposureF,

576 objects: Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None],

577 mask_plane_name: str = "INJECTED",

578 calib_flux_radius: float = 12.0,

579 draw_size_max: int = 1000,

580 inject_variance: bool = True,

581 add_noise: bool = True,

582 noise_seed: int = 0,

583 bad_mask_names: list[str] | None = None,

584 logger: Any | None = None,

585) -> tuple[list[int], list[galsim.BoundsI], list[bool], list[bool]]:

586 """Inject sources into given exposure using GalSim.

587

588 Parameters

589 ----------

590 exposure : `lsst.afw.image.ExposureF`

591 The exposure to inject synthetic sources into.

592 objects : `Generator` [`tuple`, None, None]

593 An iterator of tuples that contains (or generates) locations and object

594 surface brightness profiles to inject. The tuples should contain the

595 following elements: `lsst.geom.SpherePoint`, `lsst.geom.Point2D`,

596 `int`, `galsim.gsobject.GSObject`.

597 mask_plane_name : `str`

598 Name of the mask plane to use for the injected sources.

599 calib_flux_radius : `float`

600 Radius in pixels to use for the flux calibration. This is used to

601 produce the correct instrumental fluxes within the radius. The value

602 should match that of the field defined in slot_CalibFlux_instFlux.

603 draw_size_max : `int`

604 Maximum allowed size of the drawn object. If the object is larger than

605 this, the draw size will be clipped to this size.

606 inject_variance : `bool`

607 Whether, when injecting flux into the image plane, to inject a

608 corresponding amount of variance into the variance plane.

609 add_noise : `bool`

610 Whether to randomly vary the amount of injected image flux in each

611 pixel by an amount consistent with the amount of injected variance.

612 noise_seed : `int`

613 Seed for generating the noise of the first injected object. The seed

614 actually used increases by 1 for each subsequent object, to ensure

615 independent noise realizations.

616 bad_mask_names : `list[str]`, optional

617 List of mask plane names indicating pixels to ignore when fitting flux

618 vs variance in preparation for variance plane modification. If None,

619 then the all pixels are used.

620 logger : `lsst.utils.logging.LsstLogAdapter`, optional

621 Logger to use for logging messages.

622

623 Returns

624 -------

625 draw_sizes : `list` [`int`]

626 Draw sizes of the injected sources.

627 common_bounds : `list` [`galsim.BoundsI`]

628 Common bounds of the drawn objects.

629 fft_size_errors : `list` [`bool`]

630 Boolean flags indicating whether a GalSimFFTSizeError was raised.

631 psf_compute_errors : `list` [`bool`]

632 Boolean flags indicating whether a PSF computation error was raised.

633 """

634 exposure.mask.addMaskPlane(mask_plane_name)

635 mask_plane_core_name = mask_plane_name + "_CORE"

636 exposure.mask.addMaskPlane(mask_plane_core_name)

637 if logger:

638 logger.info(

639 "Adding %s and %s mask planes to the exposure.",

640 mask_plane_name,

641 mask_plane_core_name,

642 )

643 psf = exposure.getPsf()

644 wcs = exposure.getWcs()

645 bbox = exposure.getBBox()

646 full_bounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY)

647 galsim_image = galsim.Image(exposure.image.array, bounds=full_bounds)

648 galsim_variance = galsim.Image(exposure.variance.array, bounds=full_bounds)

649 pixel_scale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()

650

651 gain_map, is_single_CCD = get_gain_map(exposure, full_bounds, bad_mask_names)

652

653 draw_sizes: list[int] = []

654 common_bounds: list[galsim.BoundsI] = []

655 fft_size_errors: list[bool] = []

656 psf_compute_errors: list[bool] = []

657 for i, (sky_coords, pixel_coords, draw_size, object) in enumerate(objects):

658 # Instantiate default returns in case of early exit from this loop.

659 draw_sizes.append(0)

660 common_bounds.append(galsim.BoundsI())

661 fft_size_errors.append(False)

662 psf_compute_errors.append(False)

663

664 # Get spatial coordinates and WCS.

665 posd = galsim.PositionD(pixel_coords.x, pixel_coords.y)

666 posi = galsim.PositionI(pixel_coords.x // 1, pixel_coords.y // 1)

667 if logger:

668 logger.debug(f"Injecting synthetic source at {pixel_coords}.")

669 mat = wcs.linearizePixelToSky(sky_coords, arcseconds).getMatrix()

670 galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1])

671

672 # This check is here because sometimes the WCS is multivalued and

673 # objects that should not be included were being included.

674 galsim_pixel_scale = np.sqrt(galsim_wcs.pixelArea())

675 if galsim_pixel_scale < pixel_scale / 2 or galsim_pixel_scale > pixel_scale * 2:

676 continue

677

678 # Compute the PSF at the object location.

679 try:

680 psf_array = psf.computeKernelImage(pixel_coords).array

681 except InvalidParameterError:

682 # Try mapping to nearest point contained in bbox.

683 contained_point = Point2D(

684 np.clip(pixel_coords.x, bbox.minX, bbox.maxX), np.clip(pixel_coords.y, bbox.minY, bbox.maxY)

685 )

686 if pixel_coords == contained_point: # no difference, so skip immediately

687 psf_compute_errors[i] = True

688 if logger:

689 logger.debug("Cannot compute PSF for object at %s; flagging and skipping.", sky_coords)

690 continue

691 # Otherwise, try again with new point.

692 try:

693 psf_array = psf.computeKernelImage(contained_point).array

694 except InvalidParameterError:

695 psf_compute_errors[i] = True

696 if logger:

697 logger.debug("Cannot compute PSF for object at %s; flagging and skipping.", sky_coords)

698 continue

699

700 # Compute the aperture corrected PSF interpolated image.

701 aperture_correction = psf.computeApertureFlux(calib_flux_radius, psf.getAveragePosition())

702 psf_array /= aperture_correction

703 galsim_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs)

704

705 # Convolve the object with the PSF and generate draw size.

706 conv = galsim.Convolve(object, galsim_psf)

707 if draw_size == 0:

708 draw_size = conv.getGoodImageSize(galsim_wcs.minLinearScale()) # type: ignore

709 injection_draw_size = int(draw_size)

710 injection_core_size = 3

711 if draw_size_max > 0 and injection_draw_size > draw_size_max:

712 if logger:

713 logger.warning(

714 "Clipping draw size for object at %s from %d to %d pixels.",

715 sky_coords,

716 injection_draw_size,

717 draw_size_max,

718 )

719 injection_draw_size = draw_size_max

720 draw_sizes[i] = injection_draw_size

721 if injection_core_size > injection_draw_size:

722 if logger:

723 logger.debug(

724 "Clipping core size for object at %s from %d to %d pixels.",

725 sky_coords,

726 injection_core_size,

727 injection_draw_size,

728 )

729 injection_core_size = injection_draw_size

730 sub_bounds = galsim.BoundsI(posi).withBorder(injection_draw_size // 2)

731 object_common_bounds = full_bounds & sub_bounds

732 common_bounds[i] = object_common_bounds # type: ignore

733

734 # Inject the source if there is any overlap.

735 if object_common_bounds.area() > 0:

736 common_image = galsim_image[object_common_bounds]

737 common_variance = galsim_variance[object_common_bounds]

738

739 offset = posd - object_common_bounds.true_center

740 # Note, for preliminary_visit_image injection, pixel is already

741 # part of the PSF and for coadd injection, it's incorrect to

742 # include the output pixel.

743 # So for both cases, we draw using method='no_pixel'.

744 image_template = common_image.copy()

745 try:

746 image_template = conv.drawImage(

747 image_template, add_to_image=False, offset=offset, wcs=galsim_wcs, method="no_pixel"

748 )

749 except GalSimFFTSizeError as err:

750 fft_size_errors[i] = True

751 if logger:

752 logger.debug(

753 "GalSimFFTSizeError raised for object at %s; flagging and skipping.\n%s",

754 sky_coords,

755 err,

756 )

757 continue

758

759 # The amount of additional variance to inject,

760 # if inject_variance is true.

761 var_template = image_template.copy()

762 var_template *= gain_map[object_common_bounds]

763

764 if add_noise:

765 image_template, var_template = add_noise_to_galsim_image(

766 image_template,

767 var_template,

768 is_single_CCD,

769 gain_map,

770 object_common_bounds,

771 noise_seed,

772 logger,

773 )

774

775 common_image += image_template

776 if inject_variance:

777 common_variance += var_template

778

779 # Increment the seed so different noise is generated for different

780 # objects.

781 noise_seed += 1

782

783 common_box = Box2I(

784 Point2I(object_common_bounds.xmin, object_common_bounds.ymin),

785 Point2I(object_common_bounds.xmax, object_common_bounds.ymax),

786 )

787 bitvalue = exposure.mask.getPlaneBitMask(mask_plane_name)

788 exposure[common_box].mask.array |= bitvalue

789 # Add a 3 x 3 pixel mask centered on the object. The mask must be

790 # large enough to always identify the core/peak of the injected

791 # source, but small enough that it rarely overlaps real sources.

792 sub_bounds_core = galsim.BoundsI(posi).withBorder(injection_core_size // 2)

793 object_common_bounds_core = full_bounds & sub_bounds_core

794 if object_common_bounds_core.area() > 0:

795 common_box_core = Box2I(

796 Point2I(object_common_bounds_core.xmin, object_common_bounds_core.ymin),

797 Point2I(object_common_bounds_core.xmax, object_common_bounds_core.ymax),

798 )

799 bitvalue_core = exposure.mask.getPlaneBitMask(mask_plane_core_name)

800 exposure[common_box_core].mask.array |= bitvalue_core

801 else:

802 if logger:

803 logger.debug("No area overlap for object at %s; flagging and skipping.", sky_coords)

804

805 return draw_sizes, common_bounds, fft_size_errors, psf_compute_errors

Coverage for python / lsst / source / injection / inject_engine.py: 7%

298 statements