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 infer_gain_from_image(

395 exposure: ExposureF,

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

397 logger: Any | None = None,

398) -> float:

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

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

401

402 Parameters

403 ----------

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

405 The exposure to inject synthetic sources into.

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

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

408 vs. variance relationship.

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

410 Logger to use for logging messages.

411

412 Returns

413 -------

414 gain : `float`

415 The estimated average gain across all pixels.

416 """

417 image_arr = exposure.image.array

418 var_arr = exposure.variance.array

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

420

421 mask = exposure.mask

422 # Only look at mask names in the mask_plane_dict.

423 # Otherwise, getPlaneBitMask throws an exception.

424 mask_plane_dict = mask.getMaskPlaneDict()

425 if bad_mask_names is None:

426 bad_mask_names = []

427 missing = set(bad_mask_names) - set(mask_plane_dict)

428 if missing and logger:

429 logger.warning("Mask planes not found in exposure (skipping): %s", missing)

430 bad_mask_names = [name for name in bad_mask_names if name in mask_plane_dict]

431 if len(bad_mask_names) > 0:

432 bad_mask_bit_mask = mask.getPlaneBitMask(bad_mask_names)

433 good &= (mask.array.astype(int) & bad_mask_bit_mask) == 0

434

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

436 return 1.0 / fit[0]

437

438

439def get_gain_map(

440 exposure: ExposureF,

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

442 logger: Any | None = None,

443) -> galsim.Image:

444 """Retrieve gains for each pixel.

445

446 Parameters

447 ----------

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

449 The exposure to inject synthetic sources into.

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

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

452 vs. variance relationship.

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

454 Logger to use for logging messages.

455

456 Returns

457 -------

458 gain_map : `galsim.Image`

459 The gain to use in each pixel.

460 """

461 full_bbox = exposure.getBBox()

462 try:

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

464 bboxes = [amplifier.getBBox() for amplifier in amplifiers] if amplifiers else None

465 except AttributeError:

466 bboxes = None

467 if not bboxes:

468 bboxes = [full_bbox]

469

470 gains = [infer_gain_from_image(exposure.subset(bbox), bad_mask_names, logger) for bbox in bboxes]

471

472 full_bounds = galsim.BoundsI(full_bbox.minX, full_bbox.maxX, full_bbox.minY, full_bbox.maxY)

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

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

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

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

477

478 return gain_map

479

480

481def add_noise_to_galsim_image(

482 image_template: galsim.Image,

483 var_template: galsim.Image,

484 is_coadd: bool,

485 gain_map: galsim.Image,

486 object_common_bounds: galsim.BoundsI,

487 noise_seed: int | None = None,

488 logger: Any | None = None,

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

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

491 accordingly. If the image represents a coadd, add Gaussian noise.

492 Otherwise add Poisson shot noise.

493

494 Parameters

495 ----------

496 image_template : `galsim.Image`

497 The image to add noise to.

498 var_template : `galsim.Image`

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

500 is_coadd : `bool`

501 Whether the exposure is a coadd. If True, Gaussian noise is used;

502 if False, Poisson shot noise is used.

503 gain_map : `galsim.Image`

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

505 object_common_bounds : `galsim.BoundsI`

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

507 noise_seed : `int`, optional

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

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

510 Logger to use for logging messages.

511

512 Returns

513 -------

514 image_template : `galsim.Image`

515 The image with noise added.

516 var_template : `galsim.Image`

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

518 """

519 noise_template = var_template.copy()

520 negative_variance = var_template.array < 0

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

522 if np.any(negative_variance):

523 if logger:

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

525 noise_template.array[negative_variance] = 0

526 if np.any(nonfinite_variance):

527 if logger:

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

529 noise_template.array[nonfinite_variance] = 0

530

531 rng = galsim.BaseDeviate(noise_seed)

532

533 if not is_coadd:

534 # Treat noise_template, which is the expected amount of

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

536 # the number of photons collected by each pixel.

537 noise = galsim.PoissonNoise(rng)

538 noise_template.addNoise(noise)

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

540 # image flux to inject.

541 image_template = noise_template.copy()

542 image_template /= gain_map[object_common_bounds]

543 # Restore the original variance values

544 # of the bad-variance pixels.

545 bad_var = negative_variance | nonfinite_variance

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

547 var_template = noise_template

548 else:

549 noise = galsim.VariableGaussianNoise(rng, noise_template)

550 image_template.addNoise(noise)

551 var_template = image_template.copy()

552 var_template *= gain_map[object_common_bounds]

553

554 return image_template, var_template

555

556

557def inject_galsim_objects_into_exposure(

558 exposure: ExposureF,

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

560 mask_plane_name: str = "INJECTED",

561 calib_flux_radius: float = 12.0,

562 draw_size_max: int = 1000,

563 inject_variance: bool = True,

564 add_noise: bool = True,

565 noise_seed: int = 0,

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

567 logger: Any | None = None,

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

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

570

571 Parameters

572 ----------

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

574 The exposure to inject synthetic sources into.

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

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

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

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

579 `int`, `galsim.gsobject.GSObject`.

580 mask_plane_name : `str`

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

582 calib_flux_radius : `float`

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

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

585 should match that of the field defined in slot_CalibFlux_instFlux.

586 draw_size_max : `int`

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

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

589 inject_variance : `bool`

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

591 corresponding amount of variance into the variance plane.

592 add_noise : `bool`

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

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

595 noise_seed : `int`

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

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

598 independent noise realizations.

599 bad_mask_names : `list[str]`, optional

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

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

602 then the all pixels are used.

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

604 Logger to use for logging messages.

605

606 Returns

607 -------

608 draw_sizes : `list` [`int`]

609 Draw sizes of the injected sources.

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

611 Common bounds of the drawn objects.

612 fft_size_errors : `list` [`bool`]

613 Boolean flags indicating whether a GalSimFFTSizeError was raised.

614 psf_compute_errors : `list` [`bool`]

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

616 """

617 exposure.mask.addMaskPlane(mask_plane_name)

618 mask_plane_core_name = mask_plane_name + "_CORE"

619 exposure.mask.addMaskPlane(mask_plane_core_name)

620 if logger:

621 logger.info(

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

623 mask_plane_name,

624 mask_plane_core_name,

625 )

626 psf = exposure.getPsf()

627 wcs = exposure.getWcs()

628 bbox = exposure.getBBox()

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

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

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

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

633

634 gain_map = get_gain_map(exposure, bad_mask_names, logger)

635 is_coadd = exposure.info.getCoaddInputs() is not None

636

637 draw_sizes: list[int] = []

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

639 fft_size_errors: list[bool] = []

640 psf_compute_errors: list[bool] = []

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

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

643 draw_sizes.append(0)

644 common_bounds.append(galsim.BoundsI())

645 fft_size_errors.append(False)

646 psf_compute_errors.append(False)

647

648 # Get spatial coordinates and WCS.

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

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

651 if logger:

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

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

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

655

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

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

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

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

660 continue

661

662 # Compute the PSF at the object location.

663 try:

664 psf_array = psf.computeKernelImage(pixel_coords).array

665 except InvalidParameterError:

666 # Try mapping to nearest point contained in bbox.

667 contained_point = Point2D(

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

669 )

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

671 psf_compute_errors[i] = True

672 if logger:

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

674 continue

675 # Otherwise, try again with new point.

676 try:

677 psf_array = psf.computeKernelImage(contained_point).array

678 except InvalidParameterError:

679 psf_compute_errors[i] = True

680 if logger:

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

682 continue

683

684 # Compute the aperture corrected PSF interpolated image.

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

686 psf_array /= aperture_correction

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

688

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

690 conv = galsim.Convolve(object, galsim_psf)

691 if draw_size == 0:

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

693 injection_draw_size = int(draw_size)

694 injection_core_size = 3

695 if draw_size_max > 0 and injection_draw_size > draw_size_max:

696 if logger:

697 logger.warning(

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

699 sky_coords,

700 injection_draw_size,

701 draw_size_max,

702 )

703 injection_draw_size = draw_size_max

704 draw_sizes[i] = injection_draw_size

705 if injection_core_size > injection_draw_size:

706 if logger:

707 logger.debug(

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

709 sky_coords,

710 injection_core_size,

711 injection_draw_size,

712 )

713 injection_core_size = injection_draw_size

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

715 object_common_bounds = full_bounds & sub_bounds

716 common_bounds[i] = object_common_bounds # type: ignore

717

718 # Inject the source if there is any overlap.

719 if object_common_bounds.area() > 0:

720 common_image = galsim_image[object_common_bounds]

721 common_variance = galsim_variance[object_common_bounds]

722

723 offset = posd - object_common_bounds.true_center

724 # Note, for preliminary_visit_image injection, pixel is already

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

726 # include the output pixel.

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

728 image_template = common_image.copy()

729 try:

730 image_template = conv.drawImage(

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

732 )

733 except GalSimFFTSizeError as err:

734 fft_size_errors[i] = True

735 if logger:

736 logger.debug(

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

738 sky_coords,

739 err,

740 )

741 continue

742

743 # The amount of additional variance to inject,

744 # if inject_variance is true.

745 var_template = image_template.copy()

746 var_template *= gain_map[object_common_bounds]

747

748 if add_noise:

749 image_template, var_template = add_noise_to_galsim_image(

750 image_template,

751 var_template,

752 is_coadd,

753 gain_map,

754 object_common_bounds,

755 noise_seed,

756 logger,

757 )

758

759 common_image += image_template

760 if inject_variance:

761 common_variance += var_template

762

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

764 # objects.

765 noise_seed += 1

766

767 common_box = Box2I(

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

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

770 )

771 bitvalue = exposure.mask.getPlaneBitMask(mask_plane_name)

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

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

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

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

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

777 object_common_bounds_core = full_bounds & sub_bounds_core

778 if object_common_bounds_core.area() > 0:

779 common_box_core = Box2I(

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

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

782 )

783 bitvalue_core = exposure.mask.getPlaneBitMask(mask_plane_core_name)

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

785 else:

786 if logger:

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

788

789 return draw_sizes, common_bounds, fft_size_errors, psf_compute_errors

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

299 statements