Coverage for python / lsst / fgcmcal / fgcmMakeLut.py: 20%

1# See COPYRIGHT file at the top of the source tree. 

2# 

3# This file is part of fgcmcal. 

4# 

5# Developed for the LSST Data Management System. 

6# This product includes software developed by the LSST Project 

7# (https://www.lsst.org). 

8# See the COPYRIGHT file at the top-level directory of this distribution 

9# for details of code ownership. 

10# 

11# This program is free software: you can redistribute it and/or modify 

12# it under the terms of the GNU General Public License as published by 

13# the Free Software Foundation, either version 3 of the License, or 

14# (at your option) any later version. 

15# 

16# This program is distributed in the hope that it will be useful, 

17# but WITHOUT ANY WARRANTY; without even the implied warranty of 

18# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

19# GNU General Public License for more details. 

20# 

21# You should have received a copy of the GNU General Public License 

22# along with this program. If not, see <https://www.gnu.org/licenses/>. 

23"""Make a look-up-table (LUT) for FGCM calibration. 

24 

25This task computes a look-up-table for the range in expected atmosphere 

26variation and variation in instrumental throughput (as tracked by the 

27transmission_filter products). By pre-computing linearized integrals, 

28the FGCM fit is orders of magnitude faster for stars with a broad range 

29of colors and observing bands, yielding precision at the 1-2 mmag level. 

30 

31Computing a LUT requires running MODTRAN or with a pre-generated 

32atmosphere table packaged with fgcm. 

33""" 

34 

35import numpy as np 

36 

37import lsst.pex.config as pexConfig 

38import lsst.pipe.base as pipeBase 

39from lsst.pipe.base import connectionTypes 

40import lsst.afw.table as afwTable 

41import lsst.afw.cameraGeom as afwCameraGeom 

42from lsst.afw.image import TransmissionCurve 

43from .utilities import lookupStaticCalibrations, countDetectors 

44from astropy.table import Table 

45import astropy.units as units 

46 

47import fgcm 

48 

49__all__ = ['FgcmMakeLutParametersConfig', 'FgcmMakeLutConfig', 'FgcmMakeLutTask', 'SensorCorrectionTerms'] 

50 

51 

52class SensorCorrectionTerms(pexConfig.Config): 

53 refLambda = pexConfig.Field( 

54 doc="Reference wavelength for first-order correction terms.", 

55 dtype=float, 

56 optional=False, 

57 ) 

58 correctionTermDict = pexConfig.DictField( 

59 doc="Mapping of detector number to first-order correction term.", 

60 keytype=int, 

61 itemtype=float, 

62 default={}, 

63 ) 

64 

65 

66class FgcmMakeLutConnections(pipeBase.PipelineTaskConnections, 

67 dimensions=('instrument',), 

68 defaultTemplates={}): 

69 camera = connectionTypes.PrerequisiteInput( 

70 doc="Camera instrument", 

71 name="camera", 

72 storageClass="Camera", 

73 dimensions=("instrument",), 

74 isCalibration=True, 

75 ) 

76 

77 transmission_optics = connectionTypes.PrerequisiteInput( 

78 doc="Optics transmission curve information", 

79 name="transmission_optics", 

80 storageClass="TransmissionCurve", 

81 dimensions=("instrument",), 

82 isCalibration=True, 

83 deferLoad=True, 

84 ) 

85 

86 transmission_sensor = connectionTypes.PrerequisiteInput( 

87 doc="Sensor transmission curve information", 

88 name="transmission_sensor", 

89 storageClass="TransmissionCurve", 

90 dimensions=("instrument", "detector",), 

91 lookupFunction=lookupStaticCalibrations, 

92 isCalibration=True, 

93 deferLoad=True, 

94 multiple=True, 

95 ) 

96 

97 transmission_filter = connectionTypes.PrerequisiteInput( 

98 doc="Filter transmission curve information", 

99 name="transmission_filter", 

100 storageClass="TransmissionCurve", 

101 dimensions=("band", "instrument", "physical_filter",), 

102 lookupFunction=lookupStaticCalibrations, 

103 isCalibration=True, 

104 deferLoad=True, 

105 multiple=True, 

106 ) 

107 

108 transmission_filter_detector = connectionTypes.PrerequisiteInput( 

109 doc="Filter transmission curve per detector", 

110 name="transmission_filter_detector", 

111 storageClass="TransmissionCurve", 

112 dimensions=("instrument", "physical_filter", "detector",), 

113 lookupFunction=lookupStaticCalibrations, 

114 isCalibration=True, 

115 deferLoad=True, 

116 multiple=True, 

117 ) 

118 

119 fgcmLookUpTable = connectionTypes.Output( 

120 doc=("Atmosphere + instrument look-up-table for FGCM throughput and " 

121 "chromatic corrections."), 

122 name="fgcmLookUpTable", 

123 storageClass="Catalog", 

124 dimensions=("instrument",), 

125 ) 

126 

127 fgcmStandardAtmosphere = connectionTypes.Output( 

128 doc="Standard atmosphere used for FGCM calibration.", 

129 name="fgcm_standard_atmosphere", 

130 storageClass="TransmissionCurve", 

131 dimensions=("instrument",), 

132 ) 

133 

134 fgcmStandardPassbands = connectionTypes.Output( 

135 doc="Standard passbands used for FGCM calibration.", 

136 name="fgcm_standard_passband", 

137 storageClass="TransmissionCurve", 

138 dimensions=("instrument", "physical_filter"), 

139 multiple=True, 

140 ) 

141 

142 standardPassbands = connectionTypes.Output( 

143 doc="Standard passbands, in astropy table format.", 

144 name="standard_passband", 

145 storageClass="ArrowAstropy", 

146 dimensions=("instrument", "band"), 

147 multiple=True, 

148 ) 

149 

150 def __init__(self, *, config=None): 

151 if not config.doOpticsTransmission: 

152 del self.transmission_optics 

153 if not config.doSensorTransmission: 

154 del self.transmission_sensor 

155 if not config.doFilterDetectorTransmission: 

156 del self.transmission_filter_detector 

157 else: 

158 del self.transmission_filter 

159 

160 

161class FgcmMakeLutParametersConfig(pexConfig.Config): 

162 """Config for parameters if atmosphereTableName not available""" 

163 # TODO: When DM-16511 is done, it will be possible to get the 

164 # telescope elevation directly from the camera. 

165 elevation = pexConfig.Field( 

166 doc="Telescope elevation (m)", 

167 dtype=float, 

168 default=None, 

169 ) 

170 pmbRange = pexConfig.ListField( 

171 doc=("Barometric Pressure range (millibar) " 

172 "Recommended range depends on the site."), 

173 dtype=float, 

174 default=None, 

175 ) 

176 pmbSteps = pexConfig.Field( 

177 doc="Barometric Pressure number of steps", 

178 dtype=int, 

179 default=5, 

180 ) 

181 pwvRange = pexConfig.ListField( 

182 doc=("Precipitable Water Vapor range (mm) " 

183 "Recommended range depends on the site."), 

184 dtype=float, 

185 default=None, 

186 ) 

187 pwvSteps = pexConfig.Field( 

188 doc="Precipitable Water Vapor number of steps", 

189 dtype=int, 

190 default=15, 

191 ) 

192 o3Range = pexConfig.ListField( 

193 doc="Ozone range (dob)", 

194 dtype=float, 

195 default=[220.0, 310.0], 

196 ) 

197 o3Steps = pexConfig.Field( 

198 doc="Ozone number of steps", 

199 dtype=int, 

200 default=3, 

201 ) 

202 tauRange = pexConfig.ListField( 

203 doc="Aerosol Optical Depth range (unitless)", 

204 dtype=float, 

205 default=[0.002, 0.35], 

206 ) 

207 tauSteps = pexConfig.Field( 

208 doc="Aerosol Optical Depth number of steps", 

209 dtype=int, 

210 default=11, 

211 ) 

212 alphaRange = pexConfig.ListField( 

213 doc="Aerosol alpha range (unitless)", 

214 dtype=float, 

215 default=[0.0, 2.0], 

216 ) 

217 alphaSteps = pexConfig.Field( 

218 doc="Aerosol alpha number of steps", 

219 dtype=int, 

220 default=9, 

221 ) 

222 zenithRange = pexConfig.ListField( 

223 doc="Zenith angle range (degree)", 

224 dtype=float, 

225 default=[0.0, 70.0], 

226 ) 

227 zenithSteps = pexConfig.Field( 

228 doc="Zenith angle number of steps", 

229 dtype=int, 

230 default=21, 

231 ) 

232 # Note that the standard atmosphere parameters depend on the observatory 

233 # and elevation, and so these should be set on a per-camera basis. 

234 pmbStd = pexConfig.Field( 

235 doc=("Standard Atmosphere pressure (millibar); " 

236 "Recommended default depends on the site."), 

237 dtype=float, 

238 default=None, 

239 ) 

240 pwvStd = pexConfig.Field( 

241 doc=("Standard Atmosphere PWV (mm); " 

242 "Recommended default depends on the site."), 

243 dtype=float, 

244 default=None, 

245 ) 

246 o3Std = pexConfig.Field( 

247 doc="Standard Atmosphere O3 (dob)", 

248 dtype=float, 

249 default=263.0, 

250 ) 

251 tauStd = pexConfig.Field( 

252 doc="Standard Atmosphere aerosol optical depth", 

253 dtype=float, 

254 default=0.03, 

255 ) 

256 alphaStd = pexConfig.Field( 

257 doc="Standard Atmosphere aerosol alpha", 

258 dtype=float, 

259 default=1.0, 

260 ) 

261 airmassStd = pexConfig.Field( 

262 doc=("Standard Atmosphere airmass; " 

263 "Recommended default depends on the survey strategy."), 

264 dtype=float, 

265 default=None, 

266 ) 

267 lambdaNorm = pexConfig.Field( 

268 doc="Aerosol Optical Depth normalization wavelength (Angstrom)", 

269 dtype=float, 

270 default=7750.0, 

271 ) 

272 lambdaStep = pexConfig.Field( 

273 doc="Wavelength step for generating atmospheres (nm)", 

274 dtype=float, 

275 default=0.5, 

276 ) 

277 lambdaRange = pexConfig.ListField( 

278 doc="Wavelength range for LUT (Angstrom)", 

279 dtype=float, 

280 default=[3000.0, 11000.0], 

281 ) 

282 

283 

284class FgcmMakeLutConfig(pipeBase.PipelineTaskConfig, 

285 pipelineConnections=FgcmMakeLutConnections): 

286 """Config for FgcmMakeLutTask""" 

287 physicalFilters = pexConfig.ListField( 

288 doc="List of physicalFilter labels to generate look-up table.", 

289 dtype=str, 

290 default=[], 

291 ) 

292 stdPhysicalFilterOverrideMap = pexConfig.DictField( 

293 doc=("Override mapping from physical filter labels to 'standard' physical " 

294 "filter labels. The 'standard' physical filter defines the transmission " 

295 "curve that the FGCM standard bandpass will be based on. " 

296 "Any filter not listed here will be mapped to " 

297 "itself (e.g. g->g or HSC-G->HSC-G). Use this override for cross-" 

298 "filter calibration such as HSC-R->HSC-R2 and HSC-I->HSC-I2."), 

299 keytype=str, 

300 itemtype=str, 

301 default={}, 

302 ) 

303 atmosphereTableName = pexConfig.Field( 

304 doc="FGCM name or filename of precomputed atmospheres", 

305 dtype=str, 

306 default=None, 

307 optional=True, 

308 ) 

309 useScienceDetectors = pexConfig.Field( 

310 doc="Only use science detectors in LUT?", 

311 dtype=bool, 

312 default=True, 

313 ) 

314 doOpticsTransmission = pexConfig.Field( 

315 doc="Include optics transmission?", 

316 dtype=bool, 

317 default=True, 

318 ) 

319 doSensorTransmission = pexConfig.Field( 

320 doc="Include sensor transmission?", 

321 dtype=bool, 

322 default=True, 

323 ) 

324 doFilterDetectorTransmission = pexConfig.Field( 

325 doc="Use filter transmissions that are specified per-detector, rather " 

326 "than a constant or radially-dependent filter transmission?", 

327 dtype=bool, 

328 default=False, 

329 ) 

330 sensorCorrectionTermDict = pexConfig.ConfigDictField( 

331 doc="Mapping of filter name to sensor correction terms.", 

332 keytype=str, 

333 itemtype=SensorCorrectionTerms, 

334 default={}, 

335 ) 

336 parameters = pexConfig.ConfigField( 

337 doc="Atmosphere parameters (required if no atmosphereTableName)", 

338 dtype=FgcmMakeLutParametersConfig, 

339 default=None, 

340 check=None) 

341 

342 def validate(self): 

343 """ 

344 Validate the config parameters. 

345 

346 This method behaves differently from the parent validate in the case 

347 that atmosphereTableName is set. In this case, the config values 

348 for standard values, step sizes, and ranges are loaded 

349 directly from the specified atmosphereTableName. 

350 """ 

351 # check that filterNames and stdFilterNames are okay 

352 self._fields['physicalFilters'].validate(self) 

353 self._fields['stdPhysicalFilterOverrideMap'].validate(self) 

354 

355 if self.atmosphereTableName is None: 

356 # Validate the parameters 

357 self._fields['parameters'].validate(self) 

358 

359 

360class FgcmMakeLutTask(pipeBase.PipelineTask): 

361 """ 

362 Make Look-Up Table for FGCM. 

363 

364 This task computes a look-up-table for the range in expected atmosphere 

365 variation and variation in instrumental throughput (as tracked by the 

366 transmission_filter products). By pre-computing linearized integrals, 

367 the FGCM fit is orders of magnitude faster for stars with a broad range 

368 of colors and observing bands, yielding precision at the 1-2 mmag level. 

369 

370 Computing a LUT requires running MODTRAN or with a pre-generated 

371 atmosphere table packaged with fgcm. 

372 """ 

373 

374 ConfigClass = FgcmMakeLutConfig 

375 _DefaultName = "fgcmMakeLut" 

376 

377 def __init__(self, initInputs=None, **kwargs): 

378 super().__init__(**kwargs) 

379 

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

381 camera = butlerQC.get(inputRefs.camera) 

382 

383 if self.config.doOpticsTransmission: 

384 opticsHandle = butlerQC.get(inputRefs.transmission_optics) 

385 else: 

386 opticsHandle = None 

387 

388 if self.config.doSensorTransmission: 

389 sensorHandles = butlerQC.get(inputRefs.transmission_sensor) 

390 sensorHandleDict = {sensorHandle.dataId['detector']: sensorHandle for 

391 sensorHandle in sensorHandles} 

392 else: 

393 sensorHandleDict = {} 

394 

395 if self.config.doFilterDetectorTransmission: 

396 filterHandles = butlerQC.get(inputRefs.transmission_filter_detector) 

397 filterHandleDict = { 

398 (filterHandle.dataId["physical_filter"], filterHandle.dataId["detector"]): filterHandle 

399 for filterHandle in filterHandles 

400 } 

401 else: 

402 filterHandles = butlerQC.get(inputRefs.transmission_filter) 

403 filterHandleDict = {filterHandle.dataId['physical_filter']: filterHandle for 

404 filterHandle in filterHandles} 

405 

406 filterToBand = { 

407 filterHandle.dataId["physical_filter"]: filterHandle.dataId["band"] 

408 for filterHandle in filterHandles 

409 } 

410 

411 struct = self._fgcmMakeLut( 

412 camera, 

413 opticsHandle, 

414 sensorHandleDict, 

415 filterHandleDict, 

416 filterToBand, 

417 ) 

418 

419 butlerQC.put(struct.fgcmLookUpTable, outputRefs.fgcmLookUpTable) 

420 butlerQC.put(struct.fgcmStandardAtmosphere, outputRefs.fgcmStandardAtmosphere) 

421 

422 refDict = {passbandRef.dataId['physical_filter']: passbandRef for 

423 passbandRef in outputRefs.fgcmStandardPassbands} 

424 for physical_filter, passband in struct.fgcmStandardPassbands.items(): 

425 butlerQC.put(passband, refDict[physical_filter]) 

426 

427 bandRefDict = {passbandRef.dataId["band"]: passbandRef for 

428 passbandRef in outputRefs.standardPassbands} 

429 for band, passband in struct.standardPassbands.items(): 

430 butlerQC.put(passband, bandRefDict[band]) 

431 

432 def _fgcmMakeLut(self, camera, opticsHandle, sensorHandleDict, 

433 filterHandleDict, filterToBand): 

434 """ 

435 Make a FGCM Look-up Table 

436 

437 Parameters 

438 ---------- 

439 camera : `lsst.afw.cameraGeom.Camera` 

440 Camera from the butler. 

441 opticsHandle : `lsst.daf.butler.DeferredDatasetHandle` 

442 Reference to optics transmission curve. 

443 sensorHandleDict : `dict` of [`int`, `lsst.daf.butler.DeferredDatasetHandle`] 

444 Dictionary of references to sensor transmission curves. Key will 

445 be detector id. 

446 filterHandleDict : `dict` of [`str`, `lsst.daf.butler.DeferredDatasetHandle`] 

447 Dictionary of references to filter transmission curves. Key will 

448 be physical filter label or tuple of physical filter label and 

449 detector. 

450 filterToBand : `dict` [`str`: `str`] 

451 Mapping of physical filter to band name. 

452 

453 Returns 

454 ------- 

455 retStruct : `lsst.pipe.base.Struct` 

456 Output structure with keys: 

457 

458 fgcmLookUpTable : `BaseCatalog` 

459 The FGCM look-up table. 

460 fgcmStandardAtmosphere : `lsst.afw.image.TransmissionCurve` 

461 Transmission curve for the FGCM standard atmosphere. 

462 fgcmStandardPassbands : `dict` [`str`, `lsst.afw.image.TransmissionCurve`] 

463 Dictionary of fgcm standard passbands, with the key as the 

464 physical filter name. 

465 standardPassbands : `dict` [`str`, `astropy.table.Table`] 

466 Dictionary of standard passbands in astropy table format, with 

467 the key as the band name. 

468 """ 

469 # number of ccds from the length of the camera iterator 

470 nCcd = countDetectors(camera, self.config.useScienceDetectors) 

471 

472 self.log.info("Found %d ccds for look-up table" % (nCcd)) 

473 

474 # Load in optics, etc. 

475 self._loadThroughputs(camera, 

476 opticsHandle, 

477 sensorHandleDict, 

478 filterHandleDict) 

479 

480 lutConfig = self._createLutConfig(nCcd) 

481 

482 # make the lut object 

483 self.log.info("Making the LUT maker object") 

484 self.fgcmLutMaker = fgcm.FgcmLUTMaker(lutConfig) 

485 

486 # generate the throughput dictionary. 

487 

488 # these will be in Angstroms 

489 # note that lambdaStep is currently in nm, because of historical 

490 # reasons in the code. Convert to Angstroms here. 

491 throughputLambda = np.arange(self.fgcmLutMaker.lambdaRange[0], 

492 self.fgcmLutMaker.lambdaRange[1]+self.fgcmLutMaker.lambdaStep*10, 

493 self.fgcmLutMaker.lambdaStep*10.) 

494 

495 self.log.info("Built throughput lambda, %.1f-%.1f, step %.2f" % 

496 (throughputLambda[0], throughputLambda[-1], 

497 throughputLambda[1] - throughputLambda[0])) 

498 

499 throughputDict = {} 

500 for i, physicalFilter in enumerate(self.config.physicalFilters): 

501 tDict = {} 

502 tDict['LAMBDA'] = throughputLambda 

503 for ccdIndex, detector in enumerate(camera): 

504 if self.config.useScienceDetectors: 

505 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE: 

506 continue 

507 tDict[ccdIndex] = self._getThroughputDetector(detector, physicalFilter, throughputLambda) 

508 throughputDict[physicalFilter] = tDict 

509 

510 # set the throughputs 

511 self.fgcmLutMaker.setThroughputs(throughputDict) 

512 

513 # make the LUT 

514 self.log.info("Making LUT") 

515 self.fgcmLutMaker.makeLUT() 

516 

517 # and save the LUT 

518 

519 # build the index values 

520 comma = ',' 

521 physicalFilterString = comma.join(self.config.physicalFilters) 

522 stdPhysicalFilterString = comma.join(self._getStdPhysicalFilterList()) 

523 

524 atmosphereTableName = 'NoTableWasUsed' 

525 if self.config.atmosphereTableName is not None: 

526 atmosphereTableName = self.config.atmosphereTableName 

527 

528 lutSchema = self._makeLutSchema(physicalFilterString, stdPhysicalFilterString, 

529 atmosphereTableName) 

530 

531 lutCat = self._makeLutCat(lutSchema, physicalFilterString, 

532 stdPhysicalFilterString, atmosphereTableName) 

533 

534 atmStd = TransmissionCurve.makeSpatiallyConstant( 

535 throughput=self.fgcmLutMaker.atmStdTrans.astype(np.float64), 

536 wavelengths=self.fgcmLutMaker.atmLambda.astype(np.float64), 

537 throughputAtMin=self.fgcmLutMaker.atmStdTrans[0], 

538 throughputAtMax=self.fgcmLutMaker.atmStdTrans[1], 

539 ) 

540 

541 fgcmStandardPassbands = {} 

542 for i, physical_filter in enumerate(self.fgcmLutMaker.filterNames): 

543 passband = self.fgcmLutMaker.throughputs[i]['THROUGHPUT_AVG']*self.fgcmLutMaker.atmStdTrans 

544 fgcmStandardPassbands[physical_filter] = TransmissionCurve.makeSpatiallyConstant( 

545 throughput=passband.astype(np.float64), 

546 wavelengths=self.fgcmLutMaker.atmLambda.astype(np.float64), 

547 throughputAtMin=passband[0], 

548 throughputAtMax=passband[-1], 

549 ) 

550 

551 standardPassbands = {} 

552 for i, physical_filter in enumerate(self.fgcmLutMaker.filterNames): 

553 if physical_filter != self.fgcmLutMaker.stdFilterNames[i]: 

554 # This filter does not map onto one of the "standard" 

555 # passbands. E.g., for HSC we have HSC-R and HSC-R2 which 

556 # both map onto HSC-R2 which sets the "standard". 

557 continue 

558 band = filterToBand[physical_filter] 

559 passband = self.fgcmLutMaker.throughputs[i]['THROUGHPUT_AVG']*self.fgcmLutMaker.atmStdTrans 

560 passbandTable = Table( 

561 { 

562 "wavelength": (self.fgcmLutMaker.atmLambda.astype(np.float64)/10.)*units.nm, 

563 "throughput": (passband*100.)*units.percent, 

564 }, 

565 ) 

566 passbandTable["wavelength"].description = "Wavelength bin centers" 

567 standardPassbands[band] = passbandTable 

568 

569 retStruct = pipeBase.Struct( 

570 fgcmLookUpTable=lutCat, 

571 fgcmStandardAtmosphere=atmStd, 

572 fgcmStandardPassbands=fgcmStandardPassbands, 

573 standardPassbands=standardPassbands, 

574 ) 

575 

576 return retStruct 

577 

578 def _getStdPhysicalFilterList(self): 

579 """Get the standard physical filter lists from config.physicalFilters 

580 and config.stdPhysicalFilterOverrideMap 

581 

582 Returns 

583 ------- 

584 stdPhysicalFilters : `list` 

585 """ 

586 override = self.config.stdPhysicalFilterOverrideMap 

587 return [override.get(physicalFilter, physicalFilter) for 

588 physicalFilter in self.config.physicalFilters] 

589 

590 def _createLutConfig(self, nCcd): 

591 """ 

592 Create the fgcmLut config dictionary 

593 

594 Parameters 

595 ---------- 

596 nCcd: `int` 

597 Number of CCDs in the camera 

598 """ 

599 

600 # create the common stub of the lutConfig 

601 lutConfig = {} 

602 lutConfig['logger'] = self.log 

603 lutConfig['filterNames'] = self.config.physicalFilters 

604 lutConfig['stdFilterNames'] = self._getStdPhysicalFilterList() 

605 lutConfig['nCCD'] = nCcd 

606 

607 # atmosphereTable already validated if available 

608 if self.config.atmosphereTableName is not None: 

609 lutConfig['atmosphereTableName'] = self.config.atmosphereTableName 

610 else: 

611 # use the regular paramters (also validated if needed) 

612 lutConfig['elevation'] = self.config.parameters.elevation 

613 lutConfig['pmbRange'] = self.config.parameters.pmbRange 

614 lutConfig['pmbSteps'] = self.config.parameters.pmbSteps 

615 lutConfig['pwvRange'] = self.config.parameters.pwvRange 

616 lutConfig['pwvSteps'] = self.config.parameters.pwvSteps 

617 lutConfig['o3Range'] = self.config.parameters.o3Range 

618 lutConfig['o3Steps'] = self.config.parameters.o3Steps 

619 lutConfig['tauRange'] = self.config.parameters.tauRange 

620 lutConfig['tauSteps'] = self.config.parameters.tauSteps 

621 lutConfig['alphaRange'] = self.config.parameters.alphaRange 

622 lutConfig['alphaSteps'] = self.config.parameters.alphaSteps 

623 lutConfig['zenithRange'] = self.config.parameters.zenithRange 

624 lutConfig['zenithSteps'] = self.config.parameters.zenithSteps 

625 lutConfig['pmbStd'] = self.config.parameters.pmbStd 

626 lutConfig['pwvStd'] = self.config.parameters.pwvStd 

627 lutConfig['o3Std'] = self.config.parameters.o3Std 

628 lutConfig['tauStd'] = self.config.parameters.tauStd 

629 lutConfig['alphaStd'] = self.config.parameters.alphaStd 

630 lutConfig['airmassStd'] = self.config.parameters.airmassStd 

631 lutConfig['lambdaRange'] = self.config.parameters.lambdaRange 

632 lutConfig['lambdaStep'] = self.config.parameters.lambdaStep 

633 lutConfig['lambdaNorm'] = self.config.parameters.lambdaNorm 

634 

635 # Add any per-filter correction term updates if necessary. 

636 # Note that sensorCTerms is the name of the config field in fgcm. 

637 if self.config.sensorCorrectionTermDict: 

638 lutConfig['sensorCTerms'] = {} 

639 for key, value in self.config.sensorCorrectionTermDict.items(): 

640 lutConfig['sensorCTerms'][key] = ( 

641 value.refLambda, 

642 dict(value.correctionTermDict), 

643 ) 

644 

645 return lutConfig 

646 

647 def _loadThroughputs(self, camera, opticsHandle, sensorHandleDict, filterHandleDict): 

648 """Internal method to load throughput data for filters 

649 

650 Parameters 

651 ---------- 

652 camera: `lsst.afw.cameraGeom.Camera` 

653 Camera from the butler 

654 opticsHandle : `lsst.daf.butler.DeferredDatasetHandle` 

655 Reference to optics transmission curve. 

656 sensorHandleDict : `dict` of [`int`, `lsst.daf.butler.DeferredDatasetHandle`] 

657 Dictionary of references to sensor transmission curves. Key will 

658 be detector id. 

659 filterHandleDict : `dict` of [`str`, `lsst.daf.butler.DeferredDatasetHandle`] 

660 Dictionary of references to filter transmission curves. Key will 

661 be physical filter label. 

662 

663 Raises 

664 ------ 

665 ValueError : Raised if configured filter name does not match any of the 

666 available filter transmission curves. 

667 """ 

668 if self.config.doOpticsTransmission: 

669 self._opticsTransmission = opticsHandle.get() 

670 else: 

671 self._opticsTransmission = TransmissionCurve.makeSpatiallyConstant( 

672 throughput=np.ones(100), 

673 wavelengths=np.linspace( 

674 self.config.parameters.lambdaRange[0], 

675 self.config.parameters.lambdaRange[1], 

676 100, 

677 ), 

678 throughputAtMin=1.0, 

679 throughputAtMax=1.0, 

680 ) 

681 

682 self._sensorsTransmission = {} 

683 for detector in camera: 

684 if self.config.useScienceDetectors: 

685 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE: 

686 continue 

687 if self.config.doSensorTransmission: 

688 self._sensorsTransmission[detector.getId()] = sensorHandleDict[detector.getId()].get() 

689 else: 

690 self._sensorsTransmission[detector.getId()] = TransmissionCurve.makeSpatiallyConstant( 

691 throughput=np.ones(100), 

692 wavelengths=np.linspace( 

693 self.config.parameters.lambdaRange[0], 

694 self.config.parameters.lambdaRange[1], 

695 100, 

696 ), 

697 throughputAtMin=1.0, 

698 throughputAtMax=1.0, 

699 ) 

700 

701 self._filtersTransmission = {} 

702 if self.config.doFilterDetectorTransmission: 

703 for physicalFilter in self.config.physicalFilters: 

704 for detector in camera: 

705 if self.config.useScienceDetectors: 

706 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE: 

707 continue 

708 key = (physicalFilter, detector.getId()) 

709 self._filtersTransmission[key] = filterHandleDict[key].get() 

710 else: 

711 for physicalFilter in self.config.physicalFilters: 

712 self._filtersTransmission[physicalFilter] = filterHandleDict[physicalFilter].get() 

713 

714 def _getThroughputDetector(self, detector, physicalFilter, throughputLambda): 

715 """Internal method to get throughput for a detector. 

716 

717 Returns the throughput at the center of the detector for a given filter. 

718 

719 Parameters 

720 ---------- 

721 detector: `lsst.afw.cameraGeom._detector.Detector` 

722 Detector on camera 

723 physicalFilter: `str` 

724 Physical filter label 

725 throughputLambda: `np.array(dtype=np.float64)` 

726 Wavelength steps (Angstrom) 

727 

728 Returns 

729 ------- 

730 throughput: `np.array(dtype=np.float64)` 

731 Throughput (max 1.0) at throughputLambda 

732 """ 

733 

734 c = detector.getCenter(afwCameraGeom.FOCAL_PLANE) 

735 c.scale(1.0/detector.getPixelSize()[0]) # Assumes x and y pixel sizes in arcsec are the same 

736 

737 throughput = self._opticsTransmission.sampleAt(position=c, 

738 wavelengths=throughputLambda) 

739 

740 throughput *= self._sensorsTransmission[detector.getId()].sampleAt(position=c, 

741 wavelengths=throughputLambda) 

742 

743 if self.config.doFilterDetectorTransmission: 

744 throughput *= self._filtersTransmission[(physicalFilter, detector.getId())].sampleAt( 

745 position=c, 

746 wavelengths=throughputLambda, 

747 ) 

748 else: 

749 throughput *= self._filtersTransmission[physicalFilter].sampleAt( 

750 position=c, 

751 wavelengths=throughputLambda, 

752 ) 

753 

754 # Clip the throughput from 0 to 1 

755 throughput = np.clip(throughput, 0.0, 1.0) 

756 

757 return throughput 

758 

759 def _makeLutSchema(self, physicalFilterString, stdPhysicalFilterString, 

760 atmosphereTableName): 

761 """ 

762 Make the LUT schema 

763 

764 Parameters 

765 ---------- 

766 physicalFilterString: `str` 

767 Combined string of all the physicalFilters 

768 stdPhysicalFilterString: `str` 

769 Combined string of all the standard physicalFilters 

770 atmosphereTableName: `str` 

771 Name of the atmosphere table used to generate LUT 

772 

773 Returns 

774 ------- 

775 lutSchema: `afwTable.schema` 

776 """ 

777 

778 lutSchema = afwTable.Schema() 

779 

780 lutSchema.addField('tablename', type=str, doc='Atmosphere table name', 

781 size=len(atmosphereTableName)) 

782 lutSchema.addField('elevation', type=float, doc="Telescope elevation used for LUT") 

783 lutSchema.addField('physicalFilters', type=str, doc='physicalFilters in LUT', 

784 size=len(physicalFilterString)) 

785 lutSchema.addField('stdPhysicalFilters', type=str, doc='Standard physicalFilters in LUT', 

786 size=len(stdPhysicalFilterString)) 

787 lutSchema.addField('pmb', type='ArrayD', doc='Barometric Pressure', 

788 size=self.fgcmLutMaker.pmb.size) 

789 lutSchema.addField('pmbFactor', type='ArrayD', doc='PMB scaling factor', 

790 size=self.fgcmLutMaker.pmb.size) 

791 lutSchema.addField('pmbElevation', type=np.float64, doc='PMB Scaling at elevation') 

792 lutSchema.addField('pwv', type='ArrayD', doc='Preciptable Water Vapor', 

793 size=self.fgcmLutMaker.pwv.size) 

794 lutSchema.addField('o3', type='ArrayD', doc='Ozone', 

795 size=self.fgcmLutMaker.o3.size) 

796 lutSchema.addField('tau', type='ArrayD', doc='Aerosol optical depth', 

797 size=self.fgcmLutMaker.tau.size) 

798 lutSchema.addField('lambdaNorm', type=np.float64, doc='AOD wavelength') 

799 lutSchema.addField('alpha', type='ArrayD', doc='Aerosol alpha', 

800 size=self.fgcmLutMaker.alpha.size) 

801 lutSchema.addField('zenith', type='ArrayD', doc='Zenith angle', 

802 size=self.fgcmLutMaker.zenith.size) 

803 lutSchema.addField('nCcd', type=np.int32, doc='Number of CCDs') 

804 

805 # and the standard values 

806 lutSchema.addField('pmbStd', type=np.float64, doc='PMB Standard') 

807 lutSchema.addField('pwvStd', type=np.float64, doc='PWV Standard') 

808 lutSchema.addField('o3Std', type=np.float64, doc='O3 Standard') 

809 lutSchema.addField('tauStd', type=np.float64, doc='Tau Standard') 

810 lutSchema.addField('alphaStd', type=np.float64, doc='Alpha Standard') 

811 lutSchema.addField('zenithStd', type=np.float64, doc='Zenith angle Standard') 

812 lutSchema.addField('lambdaRange', type='ArrayD', doc='Wavelength range', 

813 size=2) 

814 lutSchema.addField('lambdaStep', type=np.float64, doc='Wavelength step') 

815 lutSchema.addField('lambdaStd', type='ArrayD', doc='Standard Wavelength', 

816 size=len(self.fgcmLutMaker.filterNames)) 

817 lutSchema.addField('lambdaStdFilter', type='ArrayD', doc='Standard Wavelength (raw)', 

818 size=len(self.fgcmLutMaker.filterNames)) 

819 lutSchema.addField('i0Std', type='ArrayD', doc='I0 Standard', 

820 size=len(self.fgcmLutMaker.filterNames)) 

821 lutSchema.addField('i1Std', type='ArrayD', doc='I1 Standard', 

822 size=len(self.fgcmLutMaker.filterNames)) 

823 lutSchema.addField('i10Std', type='ArrayD', doc='I10 Standard', 

824 size=len(self.fgcmLutMaker.filterNames)) 

825 lutSchema.addField('i2Std', type='ArrayD', doc='I2 Standard', 

826 size=len(self.fgcmLutMaker.filterNames)) 

827 lutSchema.addField('lambdaB', type='ArrayD', doc='Wavelength for passband (no atm)', 

828 size=len(self.fgcmLutMaker.filterNames)) 

829 lutSchema.addField('atmLambda', type='ArrayD', doc='Atmosphere wavelengths (Angstrom)',