Coverage for python / lsst / pipe / tasks / insertFakes.py: 0%

1# This file is part of pipe_tasks. 

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/>. 

21 

22""" 

23Insert fakes into deepCoadds 

24""" 

25 

26__all__ = ["InsertFakesConfig", "InsertFakesTask"] 

27 

28import galsim 

29import numpy as np 

30from astropy import units as u 

31 

32import lsst.geom as geom 

33import lsst.afw.image as afwImage 

34import lsst.afw.math as afwMath 

35import lsst.pex.config as pexConfig 

36import lsst.pipe.base as pipeBase 

37 

38from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections 

39import lsst.pipe.base.connectionTypes as cT 

40from lsst.pex.exceptions import LogicError, InvalidParameterError 

41from lsst.geom import SpherePoint, radians, Box2D, Point2D 

42 

43from deprecated.sphinx import deprecated 

44 

45 

46def _add_fake_sources(exposure, objects, calibFluxRadius=12.0, logger=None): 

47 """Add fake sources to the given exposure 

48 

49 Parameters 

50 ---------- 

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

52 The exposure into which the fake sources should be added 

53 objects : `typing.Iterator` [`tuple` ['lsst.geom.SpherePoint`, `galsim.GSObject`]] 

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

55 surface brightness profiles to inject. 

56 calibFluxRadius : `float`, optional 

57 Aperture radius (in pixels) used to define the calibration for this 

58 exposure+catalog. This is used to produce the correct instrumental fluxes 

59 within the radius. The value should match that of the field defined in 

60 slot_CalibFlux_instFlux. 

61 logger : `logging.Logger`, optional 

62 Logger. 

63 """ 

64 exposure.mask.addMaskPlane("FAKE") 

65 bitmask = exposure.mask.getPlaneBitMask("FAKE") 

66 if logger: 

67 logger.info("Adding mask plane with bitmask %s", bitmask) 

68 

69 wcs = exposure.getWcs() 

70 psf = exposure.getPsf() 

71 

72 bbox = exposure.getBBox() 

73 fullBounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY) 

74 gsImg = galsim.Image(exposure.image.array, bounds=fullBounds) 

75 

76 pixScale = wcs.getPixelScale(bbox.getCenter()).asArcseconds() 

77 

78 for spt, gsObj in objects: 

79 pt = wcs.skyToPixel(spt) 

80 posd = galsim.PositionD(pt.x, pt.y) 

81 posi = galsim.PositionI(pt.x//1, pt.y//1) 

82 if logger: 

83 logger.debug("Adding fake source at %s", pt) 

84 

85 mat = wcs.linearizePixelToSky(spt, geom.arcseconds).getMatrix() 

86 gsWCS = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) 

87 

88 # This check is here because sometimes the WCS 

89 # is multivalued and objects that should not be 

90 # were being included. 

91 gsPixScale = np.sqrt(gsWCS.pixelArea()) 

92 if gsPixScale < pixScale/2 or gsPixScale > pixScale*2: 

93 continue 

94 

95 try: 

96 psfArr = psf.computeKernelImage(pt).array 

97 apCorr = psf.computeApertureFlux(calibFluxRadius, pt) 

98 except InvalidParameterError: 

99 # Try mapping to nearest point contained in bbox. 

100 contained_pt = Point2D( 

101 np.clip(pt.x, bbox.minX, bbox.maxX), 

102 np.clip(pt.y, bbox.minY, bbox.maxY) 

103 ) 

104 if pt == contained_pt: # no difference, so skip immediately 

105 if logger: 

106 logger.info("Cannot compute Psf for object at %s; skipping", pt) 

107 continue 

108 # otherwise, try again with new point 

109 try: 

110 psfArr = psf.computeKernelImage(contained_pt).array 

111 apCorr = psf.computeApertureFlux(calibFluxRadius, contained_pt) 

112 except InvalidParameterError: 

113 if logger: 

114 logger.info("Cannot compute Psf for object at %s; skipping", pt) 

115 continue 

116 

117 psfArr /= apCorr 

118 gsPSF = galsim.InterpolatedImage(galsim.Image(psfArr), wcs=gsWCS) 

119 

120 conv = galsim.Convolve(gsObj, gsPSF) 

121 stampSize = conv.getGoodImageSize(gsWCS.minLinearScale()) 

122 subBounds = galsim.BoundsI(posi).withBorder(stampSize//2) 

123 subBounds &= fullBounds 

124 

125 if subBounds.area() > 0: 

126 subImg = gsImg[subBounds] 

127 offset = posd - subBounds.true_center 

128 # Note, for calexp injection, pixel is already part of the PSF and 

129 # for coadd injection, it's incorrect to include the output pixel. 

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

131 

132 conv.drawImage( 

133 subImg, 

134 add_to_image=True, 

135 offset=offset, 

136 wcs=gsWCS, 

137 method='no_pixel' 

138 ) 

139 

140 subBox = geom.Box2I( 

141 geom.Point2I(subBounds.xmin, subBounds.ymin), 

142 geom.Point2I(subBounds.xmax, subBounds.ymax) 

143 ) 

144 exposure[subBox].mask.array |= bitmask 

145 

146 

147def _isWCSGalsimDefault(wcs, hdr): 

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

149 or if it's just the galsim default. 

150 

151 Parameters 

152 ---------- 

153 wcs : galsim.BaseWCS 

154 Potentially default WCS. 

155 hdr : galsim.fits.FitsHeader 

156 Header as read in by galsim. 

157 

158 Returns 

159 ------- 

160 isDefault : bool 

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

162 """ 

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

164 return False 

165 if hdr.get('GS_WCS') is not None: 

166 return False 

167 if hdr.get('CTYPE1', 'LINEAR') == 'LINEAR': 

168 return not any(k in hdr for k in ['CD1_1', 'CDELT1']) 

169 for wcs_type in galsim.fitswcs.fits_wcs_types: 

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

171 try: 

172 wcs_type._readHeader(hdr) 

173 return False 

174 except Exception: 

175 pass 

176 else: 

177 return not any(k in hdr for k in ['CD1_1', 'CDELT1']) 

178 

179 

180class InsertFakesConnections(PipelineTaskConnections, 

181 defaultTemplates={"coaddName": "deep", 

182 "fakesType": "fakes_"}, 

183 dimensions=("tract", "patch", "band", "skymap")): 

184 

185 image = cT.Input( 

186 doc="Image into which fakes are to be added.", 

187 name="{coaddName}Coadd", 

188 storageClass="ExposureF", 

189 dimensions=("tract", "patch", "band", "skymap") 

190 ) 

191 

192 fakeCat = cT.Input( 

193 doc="Catalog of fake sources to draw inputs from.", 

194 name="{fakesType}fakeSourceCat", 

195 storageClass="DataFrame", 

196 dimensions=("tract", "skymap") 

197 ) 

198 

199 imageWithFakes = cT.Output( 

200 doc="Image with fake sources added.", 

201 name="{fakesType}{coaddName}Coadd", 

202 storageClass="ExposureF", 

203 dimensions=("tract", "patch", "band", "skymap") 

204 ) 

205 

206 

207@deprecated( 

208 reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.", 

209 version="v28.0", 

210 category=FutureWarning, 

211) 

212class InsertFakesConfig(PipelineTaskConfig, 

213 pipelineConnections=InsertFakesConnections): 

214 """Config for inserting fake sources 

215 """ 

216 

217 # Unchanged 

218 

219 doCleanCat = pexConfig.Field( 

220 doc="If true removes bad sources from the catalog.", 

221 dtype=bool, 

222 default=True, 

223 ) 

224 

225 fakeType = pexConfig.Field( 

226 doc="What type of fake catalog to use, snapshot (includes variability in the magnitudes calculated " 

227 "from the MJD of the image), static (no variability) or filename for a user defined fits" 

228 "catalog.", 

229 dtype=str, 

230 default="static", 

231 ) 

232 

233 calibFluxRadius = pexConfig.Field( 

234 doc="Aperture radius (in pixels) that was used to define the calibration for this image+catalog. " 

235 "This will be used to produce the correct instrumental fluxes within the radius. " 

236 "This value should match that of the field defined in slot_CalibFlux_instFlux.", 

237 dtype=float, 

238 default=12.0, 

239 ) 

240 

241 coaddName = pexConfig.Field( 

242 doc="The name of the type of coadd used", 

243 dtype=str, 

244 default="deep", 

245 ) 

246 

247 doSubSelectSources = pexConfig.Field( 

248 doc="Set to True if you wish to sub select sources to be input based on the value in the column" 

249 "set in the sourceSelectionColName config option.", 

250 dtype=bool, 

251 default=False 

252 ) 

253 

254 insertImages = pexConfig.Field( 

255 doc="Insert images directly? True or False.", 

256 dtype=bool, 

257 default=False, 

258 ) 

259 

260 insertOnlyStars = pexConfig.Field( 

261 doc="Insert only stars? True or False.", 

262 dtype=bool, 

263 default=False, 

264 ) 

265 

266 doProcessAllDataIds = pexConfig.Field( 

267 doc="If True, all input data IDs will be processed, even those containing no fake sources.", 

268 dtype=bool, 

269 default=False, 

270 ) 

271 

272 trimBuffer = pexConfig.Field( 

273 doc="Size of the pixel buffer surrounding the image. Only those fake sources with a centroid" 

274 "falling within the image+buffer region will be considered for fake source injection.", 

275 dtype=int, 

276 default=100, 

277 ) 

278 

279 sourceType = pexConfig.Field( 

280 doc="The column name for the source type used in the fake source catalog.", 

281 dtype=str, 

282 default="sourceType", 

283 ) 

284 

285 fits_alignment = pexConfig.ChoiceField( 

286 doc="How should injections from FITS files be aligned?", 

287 dtype=str, 

288 allowed={ 

289 "wcs": ( 

290 "Input image will be transformed such that the local WCS in " 

291 "the FITS header matches the local WCS in the target image. " 

292 "I.e., North, East, and angular distances in the input image " 

293 "will match North, East, and angular distances in the target " 

294 "image." 

295 ), 

296 "pixel": ( 

297 "Input image will _not_ be transformed. Up, right, and pixel " 

298 "distances in the input image will match up, right and pixel " 

299 "distances in the target image." 

300 ) 

301 }, 

302 default="pixel" 

303 ) 

304 

305 # New source catalog config variables 

306 

307 ra_col = pexConfig.Field( 

308 doc="Source catalog column name for RA (in radians).", 

309 dtype=str, 

310 default="ra", 

311 ) 

312 

313 dec_col = pexConfig.Field( 

314 doc="Source catalog column name for dec (in radians).", 

315 dtype=str, 

316 default="dec", 

317 ) 

318 

319 bulge_semimajor_col = pexConfig.Field( 

320 doc="Source catalog column name for the semimajor axis (in arcseconds) " 

321 "of the bulge half-light ellipse.", 

322 dtype=str, 

323 default="bulge_semimajor", 

324 ) 

325 

326 bulge_axis_ratio_col = pexConfig.Field( 

327 doc="Source catalog column name for the axis ratio of the bulge " 

328 "half-light ellipse.", 

329 dtype=str, 

330 default="bulge_axis_ratio", 

331 ) 

332 

333 bulge_pa_col = pexConfig.Field( 

334 doc="Source catalog column name for the position angle (measured from " 

335 "North through East in degrees) of the semimajor axis of the bulge " 

336 "half-light ellipse.", 

337 dtype=str, 

338 default="bulge_pa", 

339 ) 

340 

341 bulge_n_col = pexConfig.Field( 

342 doc="Source catalog column name for the Sersic index of the bulge.", 

343 dtype=str, 

344 default="bulge_n", 

345 ) 

346 

347 disk_semimajor_col = pexConfig.Field( 

348 doc="Source catalog column name for the semimajor axis (in arcseconds) " 

349 "of the disk half-light ellipse.", 

350 dtype=str, 

351 default="disk_semimajor", 

352 ) 

353 

354 disk_axis_ratio_col = pexConfig.Field( 

355 doc="Source catalog column name for the axis ratio of the disk " 

356 "half-light ellipse.", 

357 dtype=str, 

358 default="disk_axis_ratio", 

359 ) 

360 

361 disk_pa_col = pexConfig.Field( 

362 doc="Source catalog column name for the position angle (measured from " 

363 "North through East in degrees) of the semimajor axis of the disk " 

364 "half-light ellipse.", 

365 dtype=str, 

366 default="disk_pa", 

367 ) 

368 

369 disk_n_col = pexConfig.Field( 

370 doc="Source catalog column name for the Sersic index of the disk.", 

371 dtype=str, 

372 default="disk_n", 

373 ) 

374 

375 bulge_disk_flux_ratio_col = pexConfig.Field( 

376 doc="Source catalog column name for the bulge/disk flux ratio.", 

377 dtype=str, 

378 default="bulge_disk_flux_ratio", 

379 ) 

380 

381 mag_col = pexConfig.Field( 

382 doc="Source catalog column name template for magnitudes, in the format " 

383 "``filter name``_mag_col. E.g., if this config variable is set to " 

384 "``%s_mag``, then the i-band magnitude will be searched for in the " 

385 "``i_mag`` column of the source catalog.", 

386 dtype=str, 

387 default="%s_mag" 

388 ) 

389 

390 select_col = pexConfig.Field( 

391 doc="Source catalog column name to be used to select which sources to " 

392 "add.", 

393 dtype=str, 

394 default="select", 

395 ) 

396 

397 length_col = pexConfig.Field( 

398 doc="Source catalog column name for trail length (in pixels).", 

399 dtype=str, 

400 default="trail_length", 

401 ) 

402 

403 angle_col = pexConfig.Field( 

404 doc="Source catalog column name for trail angle (in radians).", 

405 dtype=str, 

406 default="trail_angle", 

407 ) 

408 

409 # Deprecated config variables 

410 

411 raColName = pexConfig.Field( 

412 doc="RA column name used in the fake source catalog.", 

413 dtype=str, 

414 default="raJ2000", 

415 deprecated="Use `ra_col` instead." 

416 ) 

417 

418 decColName = pexConfig.Field( 

419 doc="Dec. column name used in the fake source catalog.", 

420 dtype=str, 

421 default="decJ2000", 

422 deprecated="Use `dec_col` instead." 

423 ) 

424 

425 diskHLR = pexConfig.Field( 

426 doc="Column name for the disk half light radius used in the fake source catalog.", 

427 dtype=str, 

428 default="DiskHalfLightRadius", 

429 deprecated=( 

430 "Use `disk_semimajor_col`, `disk_axis_ratio_col`, and `disk_pa_col`" 

431 " to specify disk half-light ellipse." 

432 ) 

433 ) 

434 

435 aDisk = pexConfig.Field( 

436 doc="The column name for the semi major axis length of the disk component used in the fake source" 

437 "catalog.", 

438 dtype=str, 

439 default="a_d", 

440 deprecated=( 

441 "Use `disk_semimajor_col`, `disk_axis_ratio_col`, and `disk_pa_col`" 

442 " to specify disk half-light ellipse." 

443 ) 

444 ) 

445 

446 bDisk = pexConfig.Field( 

447 doc="The column name for the semi minor axis length of the disk component.", 

448 dtype=str, 

449 default="b_d", 

450 deprecated=( 

451 "Use `disk_semimajor_col`, `disk_axis_ratio_col`, and `disk_pa_col`" 

452 " to specify disk half-light ellipse." 

453 ) 

454 ) 

455 

456 paDisk = pexConfig.Field( 

457 doc="The column name for the PA of the disk component used in the fake source catalog.", 

458 dtype=str, 

459 default="pa_disk", 

460 deprecated=( 

461 "Use `disk_semimajor_col`, `disk_axis_ratio_col`, and `disk_pa_col`" 

462 " to specify disk half-light ellipse." 

463 ) 

464 ) 

465 

466 nDisk = pexConfig.Field( 

467 doc="The column name for the sersic index of the disk component used in the fake source catalog.", 

468 dtype=str, 

469 default="disk_n", 

470 deprecated="Use `disk_n_col` instead." 

471 ) 

472 

473 bulgeHLR = pexConfig.Field( 

474 doc="Column name for the bulge half light radius used in the fake source catalog.", 

475 dtype=str, 

476 default="BulgeHalfLightRadius", 

477 deprecated=( 

478 "Use `bulge_semimajor_col`, `bulge_axis_ratio_col`, and " 

479 "`bulge_pa_col` to specify disk half-light ellipse." 

480 ) 

481 ) 

482 

483 aBulge = pexConfig.Field( 

484 doc="The column name for the semi major axis length of the bulge component.", 

485 dtype=str, 

486 default="a_b", 

487 deprecated=( 

488 "Use `bulge_semimajor_col`, `bulge_axis_ratio_col`, and " 

489 "`bulge_pa_col` to specify disk half-light ellipse." 

490 ) 

491 ) 

492 

493 bBulge = pexConfig.Field( 

494 doc="The column name for the semi minor axis length of the bulge component used in the fake source " 

495 "catalog.", 

496 dtype=str, 

497 default="b_b", 

498 deprecated=( 

499 "Use `bulge_semimajor_col`, `bulge_axis_ratio_col`, and " 

500 "`bulge_pa_col` to specify disk half-light ellipse." 

501 ) 

502 ) 

503 

504 paBulge = pexConfig.Field( 

505 doc="The column name for the PA of the bulge component used in the fake source catalog.", 

506 dtype=str, 

507 default="pa_bulge", 

508 deprecated=( 

509 "Use `bulge_semimajor_col`, `bulge_axis_ratio_col`, and " 

510 "`bulge_pa_col` to specify disk half-light ellipse." 

511 ) 

512 ) 

513 

514 nBulge = pexConfig.Field( 

515 doc="The column name for the sersic index of the bulge component used in the fake source catalog.", 

516 dtype=str, 

517 default="bulge_n", 

518 deprecated="Use `bulge_n_col` instead." 

519 ) 

520 

521 magVar = pexConfig.Field( 

522 doc="The column name for the magnitude calculated taking variability into account. In the format " 

523 "``filter name``magVar, e.g. imagVar for the magnitude in the i band.", 

524 dtype=str, 

525 default="%smagVar", 

526 deprecated="Use `mag_col` instead." 

527 ) 

528 

529 sourceSelectionColName = pexConfig.Field( 

530 doc="The name of the column in the input fakes catalogue to be used to determine which sources to" 

531 "add, default is none and when this is used all sources are added.", 

532 dtype=str, 

533 default="templateSource", 

534 deprecated="Use `select_col` instead." 

535 ) 

536 

537 

538@deprecated( 

539 reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.", 

540 version="v28.0", 

541 category=FutureWarning, 

542) 

543class InsertFakesTask(PipelineTask): 

544 """Insert fake objects into images. 

545 

546 Add fake stars and galaxies to the given image, read in through the dataRef. Galaxy parameters are read in 

547 from the specified file and then modelled using galsim. 

548 

549 `InsertFakesTask` has five functions that make images of the fake sources and then add them to the 

550 image. 

551 

552 `addPixCoords` 

553 Use the WCS information to add the pixel coordinates of each source. 

554 `mkFakeGalsimGalaxies` 

555 Use Galsim to make fake double sersic galaxies for each set of galaxy parameters in the input file. 

556 `mkFakeStars` 

557 Use the PSF information from the image to make a fake star using the magnitude information from the 

558 input file. 

559 `cleanCat` 

560 Remove rows of the input fake catalog which have half light radius, of either the bulge or the disk, 

561 that are 0. Also removes rows that have Sersic index outside of galsim's allowed paramters. If 

562 the config option sourceSelectionColName is set then this function limits the catalog of input fakes 

563 to only those which are True in this column. 

564 `addFakeSources` 

565 Add the fake sources to the image. 

566 

567 """ 

568 

569 _DefaultName = "insertFakes" 

570 ConfigClass = InsertFakesConfig 

571 

572 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

573 inputs = butlerQC.get(inputRefs) 

574 inputs["wcs"] = inputs["image"].getWcs() 

575 inputs["photoCalib"] = inputs["image"].getPhotoCalib() 

576 

577 outputs = self.run(**inputs) 

578 butlerQC.put(outputs, outputRefs) 

579 

580 def run(self, fakeCat, image, wcs, photoCalib): 

581 """Add fake sources to an image. 

582 

583 Parameters 

584 ---------- 

585 fakeCat : `pandas.core.frame.DataFrame` 

586 The catalog of fake sources to be input 

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

588 The image into which the fake sources should be added 

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

590 WCS to use to add fake sources 

591 photoCalib : `lsst.afw.image.photoCalib.PhotoCalib` 

592 Photometric calibration to be used to calibrate the fake sources 

593 

594 Returns 

595 ------- 

596 resultStruct : `lsst.pipe.base.struct.Struct` 

597 contains : image : `lsst.afw.image.exposure.exposure.ExposureF` 

598 

599 Notes 

600 ----- 

601 Adds pixel coordinates for each source to the fakeCat and removes objects with bulge or disk half 

602 light radius = 0 (if ``config.doCleanCat = True``). 

603 

604 Adds the ``Fake`` mask plane to the image which is then set by `addFakeSources` to mark where fake 

605 sources have been added. Uses the information in the ``fakeCat`` to make fake galaxies (using galsim) 

606 and fake stars, using the PSF models from the PSF information for the image. These are then added to 

607 the image and the image with fakes included returned. 

608 

609 The galsim galaxies are made using a double sersic profile, one for the bulge and one for the disk, 

610 this is then convolved with the PSF at that point. 

611 """ 

612 # Attach overriding wcs and photoCalib to image, but retain originals 

613 # so we can reset at the end. 

614 origWcs = image.getWcs() 

615 origPhotoCalib = image.getPhotoCalib() 

616 image.setWcs(wcs) 

617 image.setPhotoCalib(photoCalib) 

618 

619 band = image.getFilter().bandLabel 

620 fakeCat = self._standardizeColumns(fakeCat, band) 

621 

622 fakeCat = self.addPixCoords(fakeCat, image) 

623 fakeCat = self.trimFakeCat(fakeCat, image) 

624 

625 if len(fakeCat) > 0: 

626 if not self.config.insertImages: 

627 if isinstance(fakeCat[self.config.sourceType].iloc[0], str): 

628 galCheckVal = "galaxy" 

629 starCheckVal = "star" 

630 trailCheckVal = "trail" 

631 elif isinstance(fakeCat[self.config.sourceType].iloc[0], bytes): 

632 galCheckVal = b"galaxy" 

633 starCheckVal = b"star" 

634 trailCheckVal = b"trail" 

635 elif isinstance(fakeCat[self.config.sourceType].iloc[0], (int, float)): 

636 galCheckVal = 1 

637 starCheckVal = 0 

638 trailCheckVal = 2 

639 else: 

640 raise TypeError( 

641 "sourceType column does not have required type, should be str, bytes or int" 

642 ) 

643 if self.config.doCleanCat: 

644 fakeCat = self.cleanCat(fakeCat, starCheckVal) 

645 

646 generator = self._generateGSObjectsFromCatalog(image, fakeCat, galCheckVal, starCheckVal, 

647 trailCheckVal) 

648 else: 

649 generator = self._generateGSObjectsFromImages(image, fakeCat) 

650 _add_fake_sources(image, generator, calibFluxRadius=self.config.calibFluxRadius, logger=self.log) 

651 elif len(fakeCat) == 0 and self.config.doProcessAllDataIds: 

652 self.log.warning("No fakes found for this dataRef; processing anyway.") 

653 image.mask.addMaskPlane("FAKE") 

654 else: 

655 raise RuntimeError("No fakes found for this dataRef.") 

656 

657 # restore original exposure WCS and photoCalib 

658 image.setWcs(origWcs) 

659 image.setPhotoCalib(origPhotoCalib) 

660 

661 resultStruct = pipeBase.Struct(imageWithFakes=image) 

662 

663 return resultStruct 

664 

665 def _standardizeColumns(self, fakeCat, band): 

666 """Use config variables to 'standardize' the expected columns and column 

667 names in the input catalog. 

668 

669 Parameters 

670 ---------- 

671 fakeCat : `pandas.core.frame.DataFrame` 

672 The catalog of fake sources to be input 

673 band : `str` 

674 Label for the current band being processed. 

675 

676 Returns 

677 ------- 

678 outCat : `pandas.core.frame.DataFrame` 

679 The standardized catalog of fake sources 

680 """ 

681 cfg = self.config 

682 replace_dict = {} 

683 

684 def add_to_replace_dict(new_name, depr_name, std_name): 

685 if new_name in fakeCat.columns: 

686 replace_dict[new_name] = std_name 

687 elif depr_name in fakeCat.columns: 

688 replace_dict[depr_name] = std_name 

689 else: 

690 raise ValueError(f"Could not determine column for {std_name}.") 

691 

692 # Prefer new config variables over deprecated config variables. 

693 # RA, dec, and mag are always required. Do these first 

694 for new_name, depr_name, std_name in [ 

695 (cfg.ra_col, cfg.raColName, 'ra'), 

696 (cfg.dec_col, cfg.decColName, 'dec'), 

697 (cfg.mag_col%band, cfg.magVar%band, 'mag') 

698 ]: 

699 add_to_replace_dict(new_name, depr_name, std_name) 

700 # Only handle bulge/disk params if not injecting images 

701 if not cfg.insertImages and not cfg.insertOnlyStars: 

702 for new_name, depr_name, std_name in [ 

703 (cfg.bulge_n_col, cfg.nBulge, 'bulge_n'), 

704 (cfg.bulge_pa_col, cfg.paBulge, 'bulge_pa'), 

705 (cfg.disk_n_col, cfg.nDisk, 'disk_n'), 

706 (cfg.disk_pa_col, cfg.paDisk, 'disk_pa'), 

707 ]: 

708 add_to_replace_dict(new_name, depr_name, std_name) 

709 

710 if cfg.doSubSelectSources: 

711 add_to_replace_dict( 

712 cfg.select_col, 

713 cfg.sourceSelectionColName, 

714 'select' 

715 ) 

716 fakeCat = fakeCat.rename(columns=replace_dict, copy=False) 

717 

718 # Handling the half-light radius and axis-ratio are trickier, since we 

719 # moved from expecting (HLR, a, b) to expecting (semimajor, axis_ratio). 

720 # Just handle these manually. 

721 if not cfg.insertImages and not cfg.insertOnlyStars: 

722 if ( 

723 cfg.bulge_semimajor_col in fakeCat.columns 

724 and cfg.bulge_axis_ratio_col in fakeCat.columns 

725 ): 

726 fakeCat = fakeCat.rename( 

727 columns={ 

728 cfg.bulge_semimajor_col: 'bulge_semimajor', 

729 cfg.bulge_axis_ratio_col: 'bulge_axis_ratio', 

730 cfg.disk_semimajor_col: 'disk_semimajor', 

731 cfg.disk_axis_ratio_col: 'disk_axis_ratio', 

732 }, 

733 copy=False 

734 ) 

735 elif ( 

736 cfg.bulgeHLR in fakeCat.columns 

737 and cfg.aBulge in fakeCat.columns 

738 and cfg.bBulge in fakeCat.columns 

739 ): 

740 fakeCat['bulge_axis_ratio'] = ( 

741 fakeCat[cfg.bBulge]/fakeCat[cfg.aBulge] 

742 ) 

743 fakeCat['bulge_semimajor'] = ( 

744 fakeCat[cfg.bulgeHLR]/np.sqrt(fakeCat['bulge_axis_ratio']) 

745 ) 

746 fakeCat['disk_axis_ratio'] = ( 

747 fakeCat[cfg.bDisk]/fakeCat[cfg.aDisk] 

748 ) 

749 fakeCat['disk_semimajor'] = ( 

750 fakeCat[cfg.diskHLR]/np.sqrt(fakeCat['disk_axis_ratio']) 

751 ) 

752 else: 

753 raise ValueError( 

754 "Could not determine columns for half-light radius and " 

755 "axis ratio." 

756 ) 

757 

758 # Process the bulge/disk flux ratio if possible. 

759 if cfg.bulge_disk_flux_ratio_col in fakeCat.columns: 

760 fakeCat = fakeCat.rename( 

761 columns={ 

762 cfg.bulge_disk_flux_ratio_col: 'bulge_disk_flux_ratio' 

763 }, 

764 copy=False 

765 ) 

766 else: 

767 fakeCat['bulge_disk_flux_ratio'] = 1.0 

768 

769 return fakeCat 

770 

771 def _generateGSObjectsFromCatalog(self, exposure, fakeCat, galCheckVal, starCheckVal, trailCheckVal): 

772 """Process catalog to generate `galsim.GSObject` s. 

773 

774 Parameters 

775 ---------- 

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

777 The exposure into which the fake sources should be added 

778 fakeCat : `pandas.core.frame.DataFrame` 

779 The catalog of fake sources to be input 

780 galCheckVal : `str`, `bytes` or `int` 

781 The value that is set in the sourceType column to specify an object is a galaxy. 

782 starCheckVal : `str`, `bytes` or `int` 

783 The value that is set in the sourceType column to specify an object is a star. 

784 trailCheckVal : `str`, `bytes` or `int` 

785 The value that is set in the sourceType column to specify an object is a star 

786 

787 Yields 

788 ------ 

789 gsObjects : `generator` 

790 A generator of tuples of `lsst.geom.SpherePoint` and `galsim.GSObject`. 

791 """ 

792 wcs = exposure.getWcs() 

793 photoCalib = exposure.getPhotoCalib() 

794 

795 self.log.info("Making %d objects for insertion", len(fakeCat)) 

796 

797 for (index, row) in fakeCat.iterrows(): 

798 ra = row['ra'] 

799 dec = row['dec'] 

800 skyCoord = SpherePoint(ra, dec, radians) 

801 xy = wcs.skyToPixel(skyCoord) 

802 

803 try: 

804 flux = photoCalib.magnitudeToInstFlux(row['mag'], xy) 

805 except LogicError: 

806 continue 

807 

808 sourceType = row[self.config.sourceType] 

809 if sourceType == galCheckVal: 

810 # GalSim convention: HLR = sqrt(a * b) = a * sqrt(b / a) 

811 bulge_gs_HLR = row['bulge_semimajor']*np.sqrt(row['bulge_axis_ratio']) 

812 bulge = galsim.Sersic(n=row['bulge_n'], half_light_radius=bulge_gs_HLR) 

813 bulge = bulge.shear(q=row['bulge_axis_ratio'], beta=((90 - row['bulge_pa'])*galsim.degrees)) 

814 

815 disk_gs_HLR = row['disk_semimajor']*np.sqrt(row['disk_axis_ratio']) 

816 disk = galsim.Sersic(n=row['disk_n'], half_light_radius=disk_gs_HLR) 

817 disk = disk.shear(q=row['disk_axis_ratio'], beta=((90 - row['disk_pa'])*galsim.degrees)) 

818 

819 gal = bulge*row['bulge_disk_flux_ratio'] + disk 

820 gal = gal.withFlux(flux) 

821 

822 yield skyCoord, gal 

823 elif sourceType == starCheckVal: 

824 star = galsim.DeltaFunction() 

825 star = star.withFlux(flux) 

826 yield skyCoord, star 

827 elif sourceType == trailCheckVal: 

828 length = row['trail_length'] 

829 angle = row['trail_angle'] 

830 

831 # Make a 'thin' box to mimic a line surface brightness profile 

832 thickness = 1e-6 # Make the box much thinner than a pixel 

833 theta = galsim.Angle(angle*galsim.radians) 

834 trail = galsim.Box(length, thickness) 

835 trail = trail.rotate(theta) 

836 trail = trail.withFlux(flux*length) 

837 

838 # Galsim objects are assumed to be in sky-coordinates. Since 

839 # we want the trail to appear as defined above in image- 

840 # coordinates, we must transform the trail here. 

841 mat = wcs.linearizePixelToSky(skyCoord, geom.arcseconds).getMatrix() 

842 trail = trail.transform(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) 

843 

844 yield skyCoord, trail 

845 else: 

846 raise TypeError(f"Unknown sourceType {sourceType}") 

847 

848 def _generateGSObjectsFromImages(self, exposure, fakeCat): 

849 """Process catalog to generate `galsim.GSObject` s. 

850 

851 Parameters 

852 ---------- 

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

854 The exposure into which the fake sources should be added 

855 fakeCat : `pandas.core.frame.DataFrame` 

856 The catalog of fake sources to be input 

857 

858 Yields 

859 ------ 

860 gsObjects : `generator` 

861 A generator of tuples of `lsst.geom.SpherePoint` and `galsim.GSObject`. 

862 """ 

863 band = exposure.getFilter().bandLabel 

864 wcs = exposure.getWcs() 

865 photoCalib = exposure.getPhotoCalib() 

866 

867 self.log.info("Processing %d fake images", len(fakeCat)) 

868 

869 for (index, row) in fakeCat.iterrows(): 

870 ra = row['ra'] 

871 dec = row['dec'] 

872 skyCoord = SpherePoint(ra, dec, radians) 

873 xy = wcs.skyToPixel(skyCoord) 

874 

875 try: 

876 flux = photoCalib.magnitudeToInstFlux(row['mag'], xy) 

877 except LogicError: 

878 continue 

879 

880 imFile = row[band+"imFilename"] 

881 try: 

882 imFile = imFile.decode("utf-8") 

883 except AttributeError: 

884 pass 

885 imFile = imFile.strip() 

886 im = galsim.fits.read(imFile, read_header=True) 

887 

888 if self.config.fits_alignment == "wcs": 

889 # galsim.fits.read will always attach a WCS to its output. If it 

890 # can't find a WCS in the FITS header, then it defaults to 

891 # scale = 1.0 arcsec / pix. So if that's the scale, then we 

892 # need to check if it was explicitly set or if it's just the 

893 # default. If it's just the default then we should raise an 

894 # exception. 

895 if _isWCSGalsimDefault(im.wcs, im.header): 

896 raise RuntimeError( 

897 f"Cannot find WCS in input FITS file {imFile}" 

898 ) 

899 elif self.config.fits_alignment == "pixel": 

900 # Here we need to set im.wcs to the local WCS at the target 

901 # position. 

902 linWcs = wcs.linearizePixelToSky(skyCoord, geom.arcseconds) 

903 mat = linWcs.getMatrix() 

904 im.wcs = galsim.JacobianWCS( 

905 mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1] 

906 ) 

907 else: 

908 raise ValueError( 

909 f"Unknown fits_alignment type {self.config.fits_alignment}"