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

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", "ElectrostaticBfMock", "DefectMock", "CrosstalkCoeffMock", 

26 "TransmissionMock", "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 brighterFatterCalibType = pexConfig.ChoiceField( 

236 dtype=str, 

237 doc="Type of calibration to generate if doGenerateData=True, otherwise ignored.", 

238 allowed={ 

239 "KERNEL": "No default suspect values; only config overrides will be used.", 

240 "ELECTROSTATIC": "Use the default from the camera model (old defaults).", 

241 }, 

242 default="KERNEL", 

243 ) 

244 

245 doDeferredCharge = pexConfig.Field( 

246 dtype=bool, 

247 default=False, 

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

249 ) 

250 doCrosstalkCoeffs = pexConfig.Field( 

251 dtype=bool, 

252 default=False, 

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

254 ) 

255 doLinearizer = pexConfig.Field( 

256 dtype=bool, 

257 default=False, 

258 doc="Return linearizer dataset.", 

259 ) 

260 doDataRef = pexConfig.Field( 

261 dtype=bool, 

262 default=False, 

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

264 ) 

265 doGenerateImage = pexConfig.Field( 

266 dtype=bool, 

267 default=False, 

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

269 ) 

270 doGenerateData = pexConfig.Field( 

271 dtype=bool, 

272 default=False, 

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

274 ) 

275 doGenerateAmpDict = pexConfig.Field( 

276 dtype=bool, 

277 default=False, 

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

279 ) 

280 

281 

282class IsrMock(pipeBase.Task): 

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

284 

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

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

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

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

289 dependent on any of the actual obs package formats. 

290 """ 

291 ConfigClass = IsrMockConfig 

292 _DefaultName = "isrMock" 

293 

294 def __init__(self, **kwargs): 

295 super().__init__(**kwargs) 

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

297 # The coefficients have units adu for IsrTask and have 

298 # units electron for IsrTaskLSST. 

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

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

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

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

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

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

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

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

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

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

309 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

310 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

311 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

312 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

313 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

314 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 

315 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]) 

316 # Generic gaussian BF kernel 

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

318 [4., 16., 26., 16., 4.], 

319 [7., 26., 41., 26., 7.], 

320 [4., 16., 26., 16., 4.], 

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

322 

323 self.aN = np.array([[1., 0.5, 0.0125, 0., 0.], 

324 [0.5, 0.010, 0., 0., 0.], 

325 [0.0125, 0., 0., 0., 0.], 

326 [0, 0., 0., 0., 0.]]) * -1e-7 

327 self.aE = np.array([[1., 0.5, 0.0125, 0., 0.], 

328 [0.5, 0.010, 0., 0., 0.], 

329 [0.0125, 0., 0., 0., 0.], 

330 [0, 0., 0., 0., 0.]]) * -1e-7 

331 self.aS = self.aN 

332 self.aW = self.aE 

333 self.aVector = (self.aN, self.aS, self.aE, self.aW) 

334 

335 def run(self): 

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

337 

338 Returns 

339 ------- 

340 image : `lsst.afw.image.Exposure` 

341 Simulated ISR image with signals added. 

342 dataProduct : 

343 Simulated ISR data products. 

344 None : 

345 Returned if no valid configuration was found. 

346 

347 Raises 

348 ------ 

349 RuntimeError 

350 Raised if both doGenerateImage and doGenerateData are specified. 

351 """ 

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

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

354 elif self.config.doGenerateImage: 

355 return self.makeImage() 

356 elif self.config.doGenerateData: 

357 return self.makeData() 

358 else: 

359 return None 

360 

361 def makeData(self): 

362 """Generate simulated ISR data. 

363 

364 Currently, only the class defined crosstalk coefficient 

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

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

367 

368 Returns 

369 ------- 

370 dataProduct : 

371 Simulated ISR data product. 

372 """ 

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

374 self.config.doDeferredCharge, 

375 self.config.doDefects, 

376 self.config.doTransmissionCurve, 

377 self.config.doCrosstalkCoeffs, 

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

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

380 elif self.config.doBrighterFatter: 

381 if self.config.brighterFatterCalibType == "KERNEL": 

382 return self.makeBfKernel() 

383 if self.config.brighterFatterCalibType == "ELECTROSTATIC": 

384 return self.makeElectrostaticBf() 

385 else: 

386 raise ValueError("%s is an unknown Brighter-Fatter calibration type." % 

387 self.config.brighterFatterCalibType) 

388 elif self.config.doDeferredCharge: 

389 return self.makeDeferredChargeCalib() 

390 elif self.config.doDefects: 

391 return self.makeDefectList() 

392 elif self.config.doTransmissionCurve: 

393 return self.makeTransmissionCurve() 

394 elif self.config.doCrosstalkCoeffs: 

395 return self.crosstalkCoeffs 

396 elif self.config.doLinearizer: 

397 return self.makeLinearizer() 

398 else: 

399 return None 

400 

401 def makeBfKernel(self): 

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

403 

404 Returns 

405 ------- 

406 kernel : `numpy.ndarray` 

407 Simulated brighter-fatter kernel. 

408 """ 

409 return self.bfKernel 

410 

411 def makeElectrostaticBf(self): 

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

413 

414 Returns 

415 ------- 

416 kernel : `numpy.ndarray` 

417 Simulated brighter-fatter kernel. 

418 """ 

419 return self.aVector 

420 

421 def makeDeferredChargeCalib(self): 

422 """Generate a CTI calibration. 

423 """ 

424 

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

426 

427 def makeDefectList(self): 

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

429 

430 Returns 

431 ------- 

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

433 Simulated defect list 

434 """ 

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

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

437 

438 def makeCrosstalkCoeff(self): 

439 """Generate the simulated crosstalk coefficients. 

440 

441 Returns 

442 ------- 

443 coeffs : `numpy.ndarray` 

444 Simulated crosstalk coefficients. 

445 """ 

446 

447 return self.crosstalkCoeffs 

448 

449 def makeTransmissionCurve(self): 

450 """Generate a simulated flat transmission curve. 

451 

452 Returns 

453 ------- 

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

455 Simulated transmission curve. 

456 """ 

457 

458 return afwImage.TransmissionCurve.makeIdentity() 

459 

460 def makeLinearity(self): 

461 """Generate a linearity dataset. 

462 

463 Returns 

464 ------- 

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

466 """ 

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

468 

469 def makeImage(self): 

470 """Generate a simulated ISR image. 

471 

472 Returns 

473 ------- 

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

475 Simulated ISR image data. 

476 

477 Notes 

478 ----- 

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

480 

481 * Generating a simulated sky with noise 

482 * Adding a single Gaussian "star" 

483 * Adding the fringe signal 

484 * Multiplying the frame by the simulated flat 

485 * Adding dark current (and noise) 

486 * Adding a bias offset (and noise) 

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

488 * Simulating crosstalk by adding a scaled version of each 

489 amplifier to each other amplifier. 

490 

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

492 three formats. 

493 

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

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

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

497 amplifier name. 

498 

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

500 implemented. 

501 

502 Note that this method generates an image in the reverse 

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

504 had a series of instrument effects added to an idealized 

505 exposure. 

506 """ 

507 exposure = self.getExposure() 

508 

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

510 bbox = None 

511 if self.config.isTrimmed is True: 

512 bbox = amp.getBBox() 

513 else: 

514 bbox = amp.getRawDataBBox() 

515 

516 ampData = exposure.image[bbox] 

517 

518 if self.config.doAddSky is True: 

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

520 

521 if self.config.doAddSource is True: 

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

523 self.config.sourceFlux, 

524 self.config.sourceX, 

525 self.config.sourceY): 

526 if idx == sourceAmp: 

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

528 

529 if self.config.doAddFringe is True: 

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

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

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

533 

534 if self.config.doAddFlat is True: 

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

536 self.amplifierAddNoise(ampData, 1.0, 0.0) 

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

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

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

540 

541 if self.config.doAddDark is True: 

542 self.amplifierAddNoise(ampData, 

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

544 np.sqrt(self.config.darkRate 

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

546 

547 if self.config.doAddCrosstalk is True: 

548 ctCalib = CrosstalkCalib() 

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

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

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

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

553 # no crosstalk mask bits are set. 

554 ctCalib.subtractCrosstalk( 

555 exposure, 

556 crosstalkCoeffs=-1*self.crosstalkCoeffs, 

557 doSubtrahendMasking=True, 

558 minPixelToMask=np.inf, 

559 ignoreVariance=True, 

560 fullAmplifier=False, 

561 ) 

562 

563 for amp in exposure.getDetector(): 

564 bbox = None 

565 if self.config.isTrimmed is True: 

566 bbox = amp.getBBox() 

567 else: 

568 bbox = amp.getRawDataBBox() 

569 

570 ampData = exposure.image[bbox] 

571 

572 if self.config.doAddBias is True: 

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

574 self.config.readNoise / self.config.gain) 

575 

576 if self.config.doAddOverscan is True: 

577 oscanBBox = amp.getRawHorizontalOverscanBBox() 

578 oscanData = exposure.image[oscanBBox] 

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

580 self.config.readNoise / self.config.gain) 

581 

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

583 1.0 * self.config.overscanScale) 

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

585 1.0 * self.config.overscanScale) 

586 

587 if self.config.doGenerateAmpDict is True: 

588 expDict = dict() 

589 for amp in exposure.getDetector(): 

590 expDict[amp.getName()] = exposure 

591 return expDict 

592 else: 

593 return exposure 

594 

595 # afw primatives to construct the image structure 

596 def getCamera(self, isForAssembly=False): 

597 """Construct a test camera object. 

598 

599 Parameters 

600 ------- 

601 isForAssembly : `bool` 

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

603 orientation (all amplifiers have LL readout 

604 corner but still contains the necessary flip 

605 and offset info needed for assembly. This is 

606 needed if isLsstLike is True. If False, return 

607 a camera with bboxes flipped and offset to the 

608 correct orientation given the readout corner. 

609 

610 Returns 

611 ------- 

612 camera : `lsst.afw.cameraGeom.camera` 

613 Test camera. 

614 """ 

615 cameraWrapper = afwTestUtils.CameraWrapper( 

616 plateScale=self.config.plateScale, 

617 radialDistortion=self.config.radialDistortion, 

618 isLsstLike=self.config.isLsstLike and isForAssembly, 

619 ) 

620 camera = cameraWrapper.camera 

621 return camera 

622 

623 def getExposure(self, isTrimmed=None): 

624 """Construct a test exposure. 

625 

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

627 unlikely header keywords that can be removed during ISR 

628 processing to exercise that code. 

629 

630 Parameters 

631 ---------- 

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

633 Override the configuration isTrimmed? 

634 

635 Returns 

636 ------- 

637 exposure : `lsst.afw.exposure.Exposure` 

638 Construct exposure containing masked image of the 

639 appropriate size. 

640 """ 

641 if isTrimmed is None: 

642 _isTrimmed = self.config.isTrimmed 

643 else: 

644 _isTrimmed = isTrimmed 

645 

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

647 detector = camera[self.config.detectorIndex] 

648 image = afwUtils.makeImageFromCcd( 

649 detector, 

650 isTrimmed=_isTrimmed, 

651 showAmpGain=False, 

652 rcMarkSize=0, 

653 binSize=1, 

654 imageFactory=afwImage.ImageF, 

655 ) 

656 

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

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

659 image.assign(0.0) 

660 

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

662 exposure = afwImage.makeExposure(maskedImage) 

663 exposure.setDetector(detector) 

664 exposure.setWcs(self.getWcs()) 

665 

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

667 exposure.getInfo().setVisitInfo(visitInfo) 

668 # Set a dummy ID. 

669 exposure.getInfo().setId(12345) 

670 

671 metadata = exposure.getMetadata() 

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

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

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

675 

676 # Add the current time 

677 now = datetime.now(timezone.utc) 

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

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

680 

681 ccd = exposure.getDetector() 

682 newCcd = ccd.rebuild() 

683 newCcd.clear() 

684 readoutMap = { 

685 'LL': ReadoutCorner.LL, 

686 'LR': ReadoutCorner.LR, 

687 'UR': ReadoutCorner.UR, 

688 'UL': ReadoutCorner.UL, 

689 } 

690 for amp in ccd: 

691 newAmp = amp.rebuild() 

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

693 newAmp.setLinearityType("Polynomial") 

694 newAmp.setGain(self.config.gain) 

695 newAmp.setSuspectLevel(25000.0) 

696 newAmp.setSaturation(32000.0) 

697 readoutCorner = amp.getReadoutCorner().name 

698 

699 # Apply flips to bbox where needed 

700 imageBBox = amp.getRawDataBBox() 

701 rawBbox = amp.getRawBBox() 

702 parallelOscanBBox = amp.getRawParallelOverscanBBox() 

703 serialOscanBBox = amp.getRawSerialOverscanBBox() 

704 prescanBBox = amp.getRawPrescanBBox() 

705 

706 if self.config.isLsstLike: 

707 # This follows cameraGeom.testUtils 

708 xoffset, yoffset = amp.getRawXYOffset() 

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

710 flipx = bool(amp.getRawFlipX()) 

711 flipy = bool(amp.getRawFlipY()) 

712 if flipx: 

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

714 rawBbox.flipLR(xExt) 

715 imageBBox.flipLR(xExt) 

716 parallelOscanBBox.flipLR(xExt) 

717 serialOscanBBox.flipLR(xExt) 

718 prescanBBox.flipLR(xExt) 

719 if flipy: 

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

721 rawBbox.flipTB(yExt) 

722 imageBBox.flipTB(yExt) 

723 parallelOscanBBox.flipTB(yExt) 

724 serialOscanBBox.flipTB(yExt) 

725 prescanBBox.flipTB(yExt) 

726 if not flipx and not flipy: 

727 readoutCorner = 'LL' 

728 elif flipx and not flipy: 

729 readoutCorner = 'LR' 

730 elif flipx and flipy: 

731 readoutCorner = 'UR' 

732 elif not flipx and flipy: 

733 readoutCorner = 'UL' 

734 rawBbox.shift(offext) 

735 imageBBox.shift(offext) 

736 parallelOscanBBox.shift(offext) 

737 serialOscanBBox.shift(offext) 

738 prescanBBox.shift(offext) 

739 newAmp.setReadoutCorner(readoutMap[readoutCorner]) 

740 newAmp.setRawBBox(rawBbox) 

741 newAmp.setRawDataBBox(imageBBox) 

742 newAmp.setRawParallelOverscanBBox(parallelOscanBBox) 

743 newAmp.setRawSerialOverscanBBox(serialOscanBBox) 

744 newAmp.setRawPrescanBBox(prescanBBox) 

745 newAmp.setRawFlipX(False) 

746 newAmp.setRawFlipY(False) 

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

748 newAmp.setRawXYOffset(no_offset) 

749 

750 newCcd.append(newAmp) 

751 

752 exposure.setDetector(newCcd.finish()) 

753 

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

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

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

757 

758 return exposure 

759 

760 def getWcs(self): 

761 """Construct a dummy WCS object. 

762 

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

764 

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

766 listed in the afwTestUtils camera definition. 

767 

768 Returns 

769 ------- 

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

771 Test WCS transform. 

772 """ 

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

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

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

776 

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

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

779 exposure coordinate. 

780 

781 Parameters 

782 ---------- 

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

784 Amplifier image to use for conversions. 

785 x : `int` 

786 X-coordinate of the point to transform. 

787 y : `int` 

788 Y-coordinate of the point to transform. 

789 

790 Returns 

791 ------- 

792 u : `int` 

793 Transformed x-coordinate. 

794 v : `int` 

795 Transformed y-coordinate. 

796 

797 Notes 

798 ----- 

799 The output is transposed intentionally here, to match the 

800 internal transpose between numpy and afw.image coordinates. 

801 """ 

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

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

804 

805 return (v, u) 

806 

807 # Simple data values. 

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

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

810 

811 This method operates in the amplifier coordinate frame. 

812 

813 Parameters 

814 ---------- 

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

816 Amplifier image to operate on. 

817 mean : `float` 

818 Mean value of the Gaussian noise. 

819 sigma : `float` 

820 Sigma of the Gaussian noise. 

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

822 Random state to use instead of self.rng. 

823 """ 

824 if rng is not None: 

825 _rng = rng 

826 else: 

827 _rng = self.rng 

828 

829 ampArr = ampData.array 

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

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

832 

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

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

835 

836 This method operates in the amplifier coordinate frame. 

837 

838 Parameters 

839 ---------- 

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

841 Amplifier image to operate on. 

842 start : `float` 

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

844 end : `float` 

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

846 """ 

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

848 ampArr = ampData.array 

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