Coverage for python/lsst/ip/isr/isrMock.py: 20%

1# This file is part of ip_isr. 

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__all__ = ["IsrMockConfig", "IsrMock", "RawMock", "TrimmedRawMock", "RawDictMock", 

23 "CalibratedRawMock", "MasterMock", 

24 "BiasMock", "DarkMock", "FlatMock", "FringeMock", "UntrimmedFringeMock", 

25 "BfKernelMock", "DefectMock", "CrosstalkCoeffMock", "TransmissionMock", 

26 "MockDataContainer", "MockFringeContainer"] 

27 

28import copy 

29import numpy as np 

30import tempfile 

31import astropy.time 

32 

33from datetime import datetime, timezone 

34 

35import lsst.geom 

36import lsst.afw.geom as afwGeom 

37import lsst.afw.image as afwImage 

38from lsstDebug import getDebugFrame 

39 

40import lsst.afw.cameraGeom.utils as afwUtils 

41import lsst.afw.cameraGeom.testUtils as afwTestUtils 

42from lsst.afw.cameraGeom import ReadoutCorner 

43import lsst.pex.config as pexConfig 

44import lsst.pipe.base as pipeBase 

45from .crosstalk import CrosstalkCalib 

46from .defects import Defects 

47 

48 

49class IsrMockConfig(pexConfig.Config): 

50 """Configuration parameters for isrMock. 

51 

52 These parameters produce generic fixed position signals from 

53 various sources, and combine them in a way that matches how those 

54 signals are combined to create real data. The camera used is the 

55 test camera defined by the afwUtils code. 

56 """ 

57 # Detector parameters. "Exposure" parameters. 

58 isLsstLike = pexConfig.Field( 

59 dtype=bool, 

60 default=False, 

61 doc="If True, products have one raw image per amplifier, otherwise, one raw image per detector.", 

62 ) 

63 plateScale = pexConfig.Field( 

64 dtype=float, 

65 default=20.0, 

66 doc="Plate scale used in constructing mock camera.", 

67 ) 

68 radialDistortion = pexConfig.Field( 

69 dtype=float, 

70 default=0.925, 

71 doc="Radial distortion term used in constructing mock camera.", 

72 ) 

73 isTrimmed = pexConfig.Field( 

74 dtype=bool, 

75 default=True, 

76 doc="If True, amplifiers have been trimmed and mosaicked to remove regions outside the data BBox.", 

77 ) 

78 detectorIndex = pexConfig.Field( 

79 dtype=int, 

80 default=20, 

81 doc="Index for the detector to use. The default value uses a standard 2x4 array of amps.", 

82 ) 

83 rngSeed = pexConfig.Field( 

84 dtype=int, 

85 default=20000913, 

86 doc="Seed for random number generator used to add noise.", 

87 ) 

88 # TODO: DM-18345 Check that mocks scale correctly when gain != 1.0 

89 gain = pexConfig.Field( 

90 dtype=float, 

91 default=1.0, 

92 doc="Gain for simulated data in electron/adu.", 

93 ) 

94 readNoise = pexConfig.Field( 

95 dtype=float, 

96 default=5.0, 

97 doc="Read noise of the detector in electron.", 

98 ) 

99 expTime = pexConfig.Field( 

100 dtype=float, 

101 default=5.0, 

102 doc="Exposure time for simulated data.", 

103 ) 

104 

105 # Signal parameters 

106 skyLevel = pexConfig.Field( 

107 dtype=float, 

108 default=1000.0, 

109 doc="Background contribution to be generated from 'the sky' in " 

110 "adu (IsrTask) or electron (IsrTaskLSST).", 

111 ) 

112 sourceFlux = pexConfig.ListField( 

113 dtype=float, 

114 default=[45000.0], 

115 doc="Peak flux level of simulated 'astronomical sources' in " 

116 "adu (IsrTask) or electron (IsrTaskLSST).", 

117 ) 

118 sourceAmp = pexConfig.ListField( 

119 dtype=int, 

120 default=[0], 

121 doc="Amplifier to place simulated 'astronomical sources'.", 

122 ) 

123 sourceX = pexConfig.ListField( 

124 dtype=float, 

125 default=[50.0], 

126 doc="Peak position (in amplifier coordinates) of simulated 'astronomical sources'.", 

127 ) 

128 sourceY = pexConfig.ListField( 

129 dtype=float, 

130 default=[25.0], 

131 doc="Peak position (in amplifier coordinates) of simulated 'astronomical sources'.", 

132 ) 

133 overscanScale = pexConfig.Field( 

134 dtype=float, 

135 default=100.0, 

136 doc="Amplitude of the ramp function to add to overscan data in " 

137 "adu (IsrTask) or electron (IsrTaskLSST)", 

138 ) 

139 biasLevel = pexConfig.Field( 

140 dtype=float, 

141 default=8000.0, 

142 doc="Background contribution to be generated from the bias offset in adu.", 

143 ) 

144 darkRate = pexConfig.Field( 

145 dtype=float, 

146 default=5.0, 

147 doc="Background level contribution (in electron/s) to be generated from dark current.", 

148 ) 

149 darkTime = pexConfig.Field( 

150 dtype=float, 

151 default=5.0, 

152 doc="Exposure time for the dark current contribution.", 

153 ) 

154 flatDrop = pexConfig.Field( 

155 dtype=float, 

156 default=0.1, 

157 doc="Fractional flux drop due to flat from center to edge of detector along x-axis.", 

158 ) 

159 fringeScale = pexConfig.ListField( 

160 dtype=float, 

161 default=[200.0], 

162 doc="Peak fluxes for the components of the fringe ripple in " 

163 "adu (IsrTask) or electron (IsrTaskLSST).", 

164 ) 

165 fringeX0 = pexConfig.ListField( 

166 dtype=float, 

167 default=[-100], 

168 doc="Center position for the fringe ripples.", 

169 ) 

170 fringeY0 = pexConfig.ListField( 

171 dtype=float, 

172 default=[-0], 

173 doc="Center position for the fringe ripples.", 

174 ) 

175 

176 # Inclusion parameters 

177 doAddSky = pexConfig.Field( 

178 dtype=bool, 

179 default=True, 

180 doc="Apply 'sky' signal to output image.", 

181 ) 

182 doAddSource = pexConfig.Field( 

183 dtype=bool, 

184 default=True, 

185 doc="Add simulated source to output image.", 

186 ) 

187 doAddCrosstalk = pexConfig.Field( 

188 dtype=bool, 

189 default=False, 

190 doc="Apply simulated crosstalk to output image. This cannot be corrected by ISR, " 

191 "as detector.hasCrosstalk()==False.", 

192 ) 

193 doAddOverscan = pexConfig.Field( 

194 dtype=bool, 

195 default=True, 

196 doc="If untrimmed, add overscan ramp to overscan and data regions.", 

197 ) 

198 doAddBias = pexConfig.Field( 

199 dtype=bool, 

200 default=True, 

201 doc="Add bias signal to data.", 

202 ) 

203 doAddDark = pexConfig.Field( 

204 dtype=bool, 

205 default=True, 

206 doc="Add dark signal to data.", 

207 ) 

208 doAddFlat = pexConfig.Field( 

209 dtype=bool, 

210 default=True, 

211 doc="Add flat signal to data.", 

212 ) 

213 doAddFringe = pexConfig.Field( 

214 dtype=bool, 

215 default=True, 

216 doc="Add fringe signal to data.", 

217 ) 

218 

219 # Datasets to create and return instead of generating an image. 

220 doTransmissionCurve = pexConfig.Field( 

221 dtype=bool, 

222 default=False, 

223 doc="Return a simulated transmission curve.", 

224 ) 

225 doDefects = pexConfig.Field( 

226 dtype=bool, 

227 default=False, 

228 doc="Return a simulated defect list.", 

229 ) 

230 doBrighterFatter = pexConfig.Field( 

231 dtype=bool, 

232 default=False, 

233 doc="Return a simulated brighter-fatter kernel.", 

234 ) 

235 doDeferredCharge = pexConfig.Field( 

236 dtype=bool, 

237 default=False, 

238 doc="Return a simulated deferred charge calibration.", 

239 ) 

240 doCrosstalkCoeffs = pexConfig.Field( 

241 dtype=bool, 

242 default=False, 

243 doc="Return the matrix of crosstalk coefficients.", 

244 ) 

245 doLinearizer = pexConfig.Field( 

246 dtype=bool, 

247 default=False, 

248 doc="Return linearizer dataset.", 

249 ) 

250 doDataRef = pexConfig.Field( 

251 dtype=bool, 

252 default=False, 

253 doc="Return a simulated gen2 butler dataRef.", 

254 ) 

255 doGenerateImage = pexConfig.Field( 

256 dtype=bool, 

257 default=False, 

258 doc="Return the generated output image if True.", 

259 ) 

260 doGenerateData = pexConfig.Field( 

261 dtype=bool, 

262 default=False, 

263 doc="Return a non-image data structure if True.", 

264 ) 

265 doGenerateAmpDict = pexConfig.Field( 

266 dtype=bool, 

267 default=False, 

268 doc="Return a dict of exposure amplifiers instead of an afwImage.Exposure.", 

269 ) 

270 

271 

272class IsrMock(pipeBase.Task): 

273 """Class to generate consistent mock images for ISR testing. 

274 

275 ISR testing currently relies on one-off fake images that do not 

276 accurately mimic the full set of detector effects. This class 

277 uses the test camera/detector/amplifier structure defined in 

278 `lsst.afw.cameraGeom.testUtils` to avoid making the test data 

279 dependent on any of the actual obs package formats. 

280 """ 

281 ConfigClass = IsrMockConfig 

282 _DefaultName = "isrMock" 

283 

284 def __init__(self, **kwargs): 

285 super().__init__(**kwargs) 

286 self.rng = np.random.RandomState(self.config.rngSeed) 

287 # The coefficients have units adu for IsrTask and have 

288 # units electron for IsrTaskLSST. 

289 self.crosstalkCoeffs = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, -1e-3, 0.0, 0.0], 

290 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

291 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

292 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

293 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

294 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

295 [1e-2, 0.0, 0.0, 2.2e-2, 0.0, 0.0, 0.0, 0.0], 

296 [1e-2, 5e-3, 5e-4, 3e-3, 4e-2, 5e-3, 5e-3, 0.0]]) 

297 if getDebugFrame(self._display, "mockCrosstalkCoeffs"): 

298 self.crosstalkCoeffs = np.array([[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

299 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

300 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

301 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

302 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

303 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

304 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

305 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]) 

306 # Generic gaussian BF kernel 

307 self.bfKernel = np.array([[1., 4., 7., 4., 1.], 

308 [4., 16., 26., 16., 4.], 

309 [7., 26., 41., 26., 7.], 

310 [4., 16., 26., 16., 4.], 

311 [1., 4., 7., 4., 1.]]) / 273.0 

312 

313 def run(self): 

314 """Generate a mock ISR product, and return it. 

315 

316 Returns 

317 ------- 

318 image : `lsst.afw.image.Exposure` 

319 Simulated ISR image with signals added. 

320 dataProduct : 

321 Simulated ISR data products. 

322 None : 

323 Returned if no valid configuration was found. 

324 

325 Raises 

326 ------ 

327 RuntimeError 

328 Raised if both doGenerateImage and doGenerateData are specified. 

329 """ 

330 if self.config.doGenerateImage and self.config.doGenerateData: 

331 raise RuntimeError("Only one of doGenerateImage and doGenerateData may be specified.") 

332 elif self.config.doGenerateImage: 

333 return self.makeImage() 

334 elif self.config.doGenerateData: 

335 return self.makeData() 

336 else: 

337 return None 

338 

339 def makeData(self): 

340 """Generate simulated ISR data. 

341 

342 Currently, only the class defined crosstalk coefficient 

343 matrix, brighter-fatter kernel, a constant unity transmission 

344 curve, or a simple single-entry defect list can be generated. 

345 

346 Returns 

347 ------- 

348 dataProduct : 

349 Simulated ISR data product. 

350 """ 

351 if sum(map(bool, [self.config.doBrighterFatter, 

352 self.config.doDeferredCharge, 

353 self.config.doDefects, 

354 self.config.doTransmissionCurve, 

355 self.config.doCrosstalkCoeffs, 

356 self.config.doLinearizer])) != 1: 

357 raise RuntimeError("Only one data product can be generated at a time.") 

358 elif self.config.doBrighterFatter: 

359 return self.makeBfKernel() 

360 elif self.config.doDeferredCharge: 

361 return self.makeDeferredChargeCalib() 

362 elif self.config.doDefects: 

363 return self.makeDefectList() 

364 elif self.config.doTransmissionCurve: 

365 return self.makeTransmissionCurve() 

366 elif self.config.doCrosstalkCoeffs: 

367 return self.crosstalkCoeffs 

368 elif self.config.doLinearizer: 

369 return self.makeLinearizer() 

370 else: 

371 return None 

372 

373 def makeBfKernel(self): 

374 """Generate a simple Gaussian brighter-fatter kernel. 

375 

376 Returns 

377 ------- 

378 kernel : `numpy.ndarray` 

379 Simulated brighter-fatter kernel. 

380 """ 

381 return self.bfKernel 

382 

383 def makeDeferredChargeCalib(self): 

384 """Generate a CTI calibration. 

385 """ 

386 

387 raise NotImplementedError("Mock deferred charge is not implemented for IsrMock.") 

388 

389 def makeDefectList(self): 

390 """Generate a simple single-entry defect list. 

391 

392 Returns 

393 ------- 

394 defectList : `lsst.meas.algorithms.Defects` 

395 Simulated defect list 

396 """ 

397 return Defects([lsst.geom.Box2I(lsst.geom.Point2I(0, 0), 

398 lsst.geom.Extent2I(40, 50))]) 

399 

400 def makeCrosstalkCoeff(self): 

401 """Generate the simulated crosstalk coefficients. 

402 

403 Returns 

404 ------- 

405 coeffs : `numpy.ndarray` 

406 Simulated crosstalk coefficients. 

407 """ 

408 

409 return self.crosstalkCoeffs 

410 

411 def makeTransmissionCurve(self): 

412 """Generate a simulated flat transmission curve. 

413 

414 Returns 

415 ------- 

416 transmission : `lsst.afw.image.TransmissionCurve` 

417 Simulated transmission curve. 

418 """ 

419 

420 return afwImage.TransmissionCurve.makeIdentity() 

421 

422 def makeLinearity(self): 

423 """Generate a linearity dataset. 

424 

425 Returns 

426 ------- 

427 linearizer : `lsst.ip.isr.Linearizer` 

428 """ 

429 raise NotImplementedError("Linearizer not implemented for isrMock.") 

430 

431 def makeImage(self): 

432 """Generate a simulated ISR image. 

433 

434 Returns 

435 ------- 

436 exposure : `lsst.afw.image.Exposure` or `dict` 

437 Simulated ISR image data. 

438 

439 Notes 

440 ----- 

441 This method currently constructs a "raw" data image by: 

442 

443 * Generating a simulated sky with noise 

444 * Adding a single Gaussian "star" 

445 * Adding the fringe signal 

446 * Multiplying the frame by the simulated flat 

447 * Adding dark current (and noise) 

448 * Adding a bias offset (and noise) 

449 * Adding an overscan gradient parallel to the pixel y-axis 

450 * Simulating crosstalk by adding a scaled version of each 

451 amplifier to each other amplifier. 

452 

453 The exposure with image data constructed this way is in one of 

454 three formats. 

455 

456 * A single image, with overscan and prescan regions retained 

457 * A single image, with overscan and prescan regions trimmed 

458 * A `dict`, containing the amplifer data indexed by the 

459 amplifier name. 

460 

461 The nonlinearity, CTE, and brighter fatter are currently not 

462 implemented. 

463 

464 Note that this method generates an image in the reverse 

465 direction as the ISR processing, as the output image here has 

466 had a series of instrument effects added to an idealized 

467 exposure. 

468 """ 

469 exposure = self.getExposure() 

470 

471 for idx, amp in enumerate(exposure.getDetector()): 

472 bbox = None 

473 if self.config.isTrimmed is True: 

474 bbox = amp.getBBox() 

475 else: 

476 bbox = amp.getRawDataBBox() 

477 

478 ampData = exposure.image[bbox] 

479 

480 if self.config.doAddSky is True: 

481 self.amplifierAddNoise(ampData, self.config.skyLevel, np.sqrt(self.config.skyLevel)) 

482 

483 if self.config.doAddSource is True: 

484 for sourceAmp, sourceFlux, sourceX, sourceY in zip(self.config.sourceAmp, 

485 self.config.sourceFlux, 

486 self.config.sourceX, 

487 self.config.sourceY): 

488 if idx == sourceAmp: 

489 self.amplifierAddSource(ampData, sourceFlux, sourceX, sourceY) 

490 

491 if self.config.doAddFringe is True: 

492 self.amplifierAddFringe(amp, ampData, np.array(self.config.fringeScale), 

493 x0=np.array(self.config.fringeX0), 

494 y0=np.array(self.config.fringeY0)) 

495 

496 if self.config.doAddFlat is True: 

497 if ampData.getArray().sum() == 0.0: 

498 self.amplifierAddNoise(ampData, 1.0, 0.0) 

499 u0 = exposure.getDimensions().getX() 

500 v0 = exposure.getDimensions().getY() 

501 self.amplifierMultiplyFlat(amp, ampData, self.config.flatDrop, u0=u0, v0=v0) 

502 

503 if self.config.doAddDark is True: 

504 self.amplifierAddNoise(ampData, 

505 self.config.darkRate * self.config.darkTime / self.config.gain, 

506 np.sqrt(self.config.darkRate 

507 * self.config.darkTime / self.config.gain)) 

508 

509 if self.config.doAddCrosstalk is True: 

510 ctCalib = CrosstalkCalib() 

511 # We use the regular subtractCrosstalk code but with a negative 

512 # sign on the crosstalk coefficients so it adds instead of 

513 # subtracts. We only apply the signal plane (ignoreVariance, 

514 # subtrahendMasking) with a very large pixel to mask to ensure 

515 # no crosstalk mask bits are set. 

516 ctCalib.subtractCrosstalk( 

517 exposure, 

518 crosstalkCoeffs=-1*self.crosstalkCoeffs, 

519 doSubtrahendMasking=True, 

520 minPixelToMask=np.inf, 

521 ignoreVariance=True, 

522 fullAmplifier=False, 

523 ) 

524 

525 for amp in exposure.getDetector(): 

526 bbox = None 

527 if self.config.isTrimmed is True: 

528 bbox = amp.getBBox() 

529 else: 

530 bbox = amp.getRawDataBBox() 

531 

532 ampData = exposure.image[bbox] 

533 

534 if self.config.doAddBias is True: 

535 self.amplifierAddNoise(ampData, self.config.biasLevel, 

536 self.config.readNoise / self.config.gain) 

537 

538 if self.config.doAddOverscan is True: 

539 oscanBBox = amp.getRawHorizontalOverscanBBox() 

540 oscanData = exposure.image[oscanBBox] 

541 self.amplifierAddNoise(oscanData, self.config.biasLevel, 

542 self.config.readNoise / self.config.gain) 

543 

544 self.amplifierAddYGradient(ampData, -1.0 * self.config.overscanScale, 

545 1.0 * self.config.overscanScale) 

546 self.amplifierAddYGradient(oscanData, -1.0 * self.config.overscanScale, 

547 1.0 * self.config.overscanScale) 

548 

549 if self.config.doGenerateAmpDict is True: 

550 expDict = dict() 

551 for amp in exposure.getDetector(): 

552 expDict[amp.getName()] = exposure 

553 return expDict 

554 else: 

555 return exposure 

556 

557 # afw primatives to construct the image structure 

558 def getCamera(self, isForAssembly=False): 

559 """Construct a test camera object. 

560 

561 Parameters 

562 ------- 

563 isForAssembly : `bool` 

564 If True, construct a camera with "super raw" 

565 orientation (all amplifiers have LL readout 

566 corner but still contains the necessary flip 

567 and offset info needed for assembly. This is 

568 needed if isLsstLike is True. If False, return 

569 a camera with bboxes flipped and offset to the 

570 correct orientation given the readout corner. 

571 

572 Returns 

573 ------- 

574 camera : `lsst.afw.cameraGeom.camera` 

575 Test camera. 

576 """ 

577 cameraWrapper = afwTestUtils.CameraWrapper( 

578 plateScale=self.config.plateScale, 

579 radialDistortion=self.config.radialDistortion, 

580 isLsstLike=self.config.isLsstLike and isForAssembly, 

581 ) 

582 camera = cameraWrapper.camera 

583 return camera 

584 

585 def getExposure(self, isTrimmed=None): 

586 """Construct a test exposure. 

587 

588 The test exposure has a simple WCS set, as well as a list of 

589 unlikely header keywords that can be removed during ISR 

590 processing to exercise that code. 

591 

592 Parameters 

593 ---------- 

594 isTrimmed : `bool` or `None`, optional 

595 Override the configuration isTrimmed? 

596 

597 Returns 

598 ------- 

599 exposure : `lsst.afw.exposure.Exposure` 

600 Construct exposure containing masked image of the 

601 appropriate size. 

602 """ 

603 if isTrimmed is None: 

604 _isTrimmed = self.config.isTrimmed 

605 else: 

606 _isTrimmed = isTrimmed 

607 

608 camera = self.getCamera(isForAssembly=self.config.isLsstLike) 

609 detector = camera[self.config.detectorIndex] 

610 image = afwUtils.makeImageFromCcd( 

611 detector, 

612 isTrimmed=_isTrimmed, 

613 showAmpGain=False, 

614 rcMarkSize=0, 

615 binSize=1, 

616 imageFactory=afwImage.ImageF, 

617 ) 

618 

619 var = afwImage.ImageF(image.getDimensions()) 

620 mask = afwImage.Mask(image.getDimensions()) 

621 image.assign(0.0) 

622 

623 maskedImage = afwImage.makeMaskedImage(image, mask, var) 

624 exposure = afwImage.makeExposure(maskedImage) 

625 exposure.setDetector(detector) 

626 exposure.setWcs(self.getWcs()) 

627 

628 visitInfo = afwImage.VisitInfo(exposureTime=self.config.expTime, darkTime=self.config.darkTime) 

629 exposure.getInfo().setVisitInfo(visitInfo) 

630 # Set a dummy ID. 

631 exposure.getInfo().setId(12345) 

632 

633 metadata = exposure.getMetadata() 

634 metadata.add("SHEEP", 7.3, "number of sheep on farm") 

635 metadata.add("MONKEYS", 155, "monkeys per tree") 

636 metadata.add("VAMPIRES", 4, "How scary are vampires.") 

637 

638 # Add the current time 

639 now = datetime.now(timezone.utc) 

640 currentMjd = astropy.time.Time(now, format='datetime', scale='utc').mjd 

641 metadata.add("MJD", currentMjd, "Modified Julian Date that the file was written") 

642 

643 ccd = exposure.getDetector() 

644 newCcd = ccd.rebuild() 

645 newCcd.clear() 

646 readoutMap = { 

647 'LL': ReadoutCorner.LL, 

648 'LR': ReadoutCorner.LR, 

649 'UR': ReadoutCorner.UR, 

650 'UL': ReadoutCorner.UL, 

651 } 

652 for amp in ccd: 

653 newAmp = amp.rebuild() 

654 newAmp.setLinearityCoeffs((0., 1., 0., 0.)) 

655 newAmp.setLinearityType("Polynomial") 

656 newAmp.setGain(self.config.gain) 

657 newAmp.setSuspectLevel(25000.0) 

658 newAmp.setSaturation(32000.0) 

659 readoutCorner = amp.getReadoutCorner().name 

660 

661 # Apply flips to bbox where needed 

662 imageBBox = amp.getRawDataBBox() 

663 rawBbox = amp.getRawBBox() 

664 parallelOscanBBox = amp.getRawParallelOverscanBBox() 

665 serialOscanBBox = amp.getRawSerialOverscanBBox() 

666 prescanBBox = amp.getRawPrescanBBox() 

667 

668 if self.config.isLsstLike: 

669 # This follows cameraGeom.testUtils 

670 xoffset, yoffset = amp.getRawXYOffset() 

671 offext = lsst.geom.Extent2I(xoffset, yoffset) 

672 flipx = bool(amp.getRawFlipX()) 

673 flipy = bool(amp.getRawFlipY()) 

674 if flipx: 

675 xExt = rawBbox.getDimensions().getX() 

676 rawBbox.flipLR(xExt) 

677 imageBBox.flipLR(xExt) 

678 parallelOscanBBox.flipLR(xExt) 

679 serialOscanBBox.flipLR(xExt) 

680 prescanBBox.flipLR(xExt) 

681 if flipy: 

682 yExt = rawBbox.getDimensions().getY() 

683 rawBbox.flipTB(yExt) 

684 imageBBox.flipTB(yExt) 

685 parallelOscanBBox.flipTB(yExt) 

686 serialOscanBBox.flipTB(yExt) 

687 prescanBBox.flipTB(yExt) 

688 if not flipx and not flipy: 

689 readoutCorner = 'LL' 

690 elif flipx and not flipy: 

691 readoutCorner = 'LR' 

692 elif flipx and flipy: 

693 readoutCorner = 'UR' 

694 elif not flipx and flipy: 

695 readoutCorner = 'UL' 

696 rawBbox.shift(offext) 

697 imageBBox.shift(offext) 

698 parallelOscanBBox.shift(offext) 

699 serialOscanBBox.shift(offext) 

700 prescanBBox.shift(offext) 

701 newAmp.setReadoutCorner(readoutMap[readoutCorner]) 

702 newAmp.setRawBBox(rawBbox) 

703 newAmp.setRawDataBBox(imageBBox) 

704 newAmp.setRawParallelOverscanBBox(parallelOscanBBox) 

705 newAmp.setRawSerialOverscanBBox(serialOscanBBox) 

706 newAmp.setRawPrescanBBox(prescanBBox) 

707 newAmp.setRawFlipX(False) 

708 newAmp.setRawFlipY(False) 

709 no_offset = lsst.geom.Extent2I(0, 0) 

710 newAmp.setRawXYOffset(no_offset) 

711 

712 newCcd.append(newAmp) 

713 

714 exposure.setDetector(newCcd.finish()) 

715 

716 exposure.image.array[:] = np.zeros(exposure.getImage().getDimensions()).transpose() 

717 exposure.mask.array[:] = np.zeros(exposure.getMask().getDimensions()).transpose() 

718 exposure.variance.array[:] = np.zeros(exposure.getVariance().getDimensions()).transpose() 

719 

720 return exposure 

721 

722 def getWcs(self): 

723 """Construct a dummy WCS object. 

724 

725 Taken from the deprecated ip_isr/examples/exampleUtils.py. 

726 

727 This is not guaranteed, given the distortion and pixel scale 

728 listed in the afwTestUtils camera definition. 

729 

730 Returns 

731 ------- 

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

733 Test WCS transform. 

734 """ 

735 return afwGeom.makeSkyWcs(crpix=lsst.geom.Point2D(0.0, 100.0), 

736 crval=lsst.geom.SpherePoint(45.0, 25.0, lsst.geom.degrees), 

737 cdMatrix=afwGeom.makeCdMatrix(scale=1.0*lsst.geom.degrees)) 

738 

739 def localCoordToExpCoord(self, ampData, x, y): 

740 """Convert between a local amplifier coordinate and the full 

741 exposure coordinate. 

742 

743 Parameters 

744 ---------- 

745 ampData : `lsst.afw.image.ImageF` 

746 Amplifier image to use for conversions. 

747 x : `int` 

748 X-coordinate of the point to transform. 

749 y : `int` 

750 Y-coordinate of the point to transform. 

751 

752 Returns 

753 ------- 

754 u : `int` 

755 Transformed x-coordinate. 

756 v : `int` 

757 Transformed y-coordinate. 

758 

759 Notes 

760 ----- 

761 The output is transposed intentionally here, to match the 

762 internal transpose between numpy and afw.image coordinates. 

763 """ 

764 u = x + ampData.getBBox().getBeginX() 

765 v = y + ampData.getBBox().getBeginY() 

766 

767 return (v, u) 

768 

769 # Simple data values. 

770 def amplifierAddNoise(self, ampData, mean, sigma, rng=None): 

771 """Add Gaussian noise to an amplifier's image data. 

772 

773 This method operates in the amplifier coordinate frame. 

774 

775 Parameters 

776 ---------- 

777 ampData : `lsst.afw.image.ImageF` 

778 Amplifier image to operate on. 

779 mean : `float` 

780 Mean value of the Gaussian noise. 

781 sigma : `float` 

782 Sigma of the Gaussian noise. 

783 rng : `np.random.RandomState`, optional 

784 Random state to use instead of self.rng. 

785 """ 

786 if rng is not None: 

787 _rng = rng 

788 else: 

789 _rng = self.rng 

790 

791 ampArr = ampData.array 

792 ampArr[:] = ampArr[:] + _rng.normal(mean, sigma, 

793 size=ampData.getDimensions()).transpose() 

794 

795 def amplifierAddYGradient(self, ampData, start, end): 

796 """Add a y-axis linear gradient to an amplifier's image data. 

797 

798 This method operates in the amplifier coordinate frame. 

799 

800 Parameters 

801 ---------- 

802 ampData : `lsst.afw.image.ImageF` 

803 Amplifier image to operate on. 

804 start : `float` 

805 Start value of the gradient (at y=0). 

806 end : `float` 

807 End value of the gradient (at y=ymax). 

808 """ 

809 nPixY = ampData.getDimensions().getY() 

810 ampArr = ampData.array 

811 ampArr[:] = ampArr[:] + (np.interp(range(nPixY), (0, nPixY - 1), (start, end)).reshape(nPixY, 1) 

812 + np.zeros(ampData.getDimensions()).transpose()) 

813 

814 def amplifierAddSource(self, ampData, scale, x0, y0): 

815 """Add a single Gaussian source to an amplifier. 

816 

817 This method operates in the amplifier coordinate frame. 

818 

819 Parameters 

820 ---------- 

821 ampData : `lsst.afw.image.ImageF` 

822 Amplifier image to operate on. 

823 scale : `float` 

824 Peak flux of the source to add. 

825 x0 : `float` 

826 X-coordinate of the source peak. 

827 y0 : `float` 

828 Y-coordinate of the source peak. 

829 """ 

830 for x in range(0, ampData.getDimensions().getX()): 

831 for y in range(0, ampData.getDimensions().getY()): 

832 ampData.array[y][x] = (ampData.array[y][x] 

833 + scale * np.exp(-0.5 * ((x - x0)**2 + (y - y0)**2) / 3.0**2)) 

834 

835 # Functional form data values. 

836 def amplifierAddFringe(self, amp, ampData, scale, x0=100, y0=0): 

837 """Add a fringe-like ripple pattern to an amplifier's image data. 

838 

839 Parameters 

840 ---------- 

841 amp : `~lsst.afw.ampInfo.AmpInfoRecord` 

842 Amplifier to operate on. Needed for amp<->exp coordinate 

843 transforms. 

844 ampData : `lsst.afw.image.ImageF` 

845 Amplifier image to operate on. 

846 scale : `numpy.array` or `float` 

847 Peak intensity scaling for the ripple. 

848 x0 : `numpy.array` or `float`, optional 

849 Fringe center 

850 y0 : `numpy.array` or `float`, optional 

851 Fringe center 

852 

853 Notes 

854 ----- 

855 This uses an offset sinc function to generate a ripple 

856 pattern. True fringes have much finer structure, but this 

857 pattern should be visually identifiable. The (x, y)