Coverage for python / lsst / fgcmcal / fgcmOutputProducts.py: 14%

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 the final fgcmcal output products. 

24 

25This task takes the final output from fgcmFitCycle and produces the following 

26outputs for use in the DM stack: the FGCM standard stars in a reference 

27catalog format; the model atmospheres in "transmission_atmosphere_fgcm" 

28format; and the zeropoints in "fgcm_photoCalib" format. Optionally, the 

29task can transfer the 'absolute' calibration from a reference catalog 

30to put the fgcm standard stars in units of Jansky. This is accomplished 

31by matching stars in a sample of healpix pixels, and applying the median 

32offset per band. 

33""" 

34import copy 

35 

36import numpy as np 

37import hpgeom as hpg 

38from astropy import units 

39from astropy.table import Table 

40import esutil 

41 

42import lsst.daf.base as dafBase 

43import lsst.pex.config as pexConfig 

44import lsst.pipe.base as pipeBase 

45from lsst.pipe.base import connectionTypes 

46from lsst.afw.image import TransmissionCurve 

47from lsst.meas.algorithms import ReferenceObjectLoader, LoadReferenceObjectsConfig 

48from lsst.pipe.tasks.photoCal import PhotoCalTask 

49import lsst.geom 

50import lsst.afw.image as afwImage 

51import lsst.afw.math as afwMath 

52import lsst.afw.table as afwTable 

53from lsst.skymap import BaseSkyMap 

54 

55from .utilities import computeApproxPixelAreaFields 

56from .utilities import FGCM_ILLEGAL_VALUE 

57 

58import fgcm 

59 

60__all__ = ['FgcmOutputProductsConfig', 'FgcmOutputProductsTask'] 

61 

62 

63class FgcmOutputProductsConnections(pipeBase.PipelineTaskConnections, 

64 dimensions=("instrument",), 

65 defaultTemplates={"cycleNumber": "0"}): 

66 camera = connectionTypes.PrerequisiteInput( 

67 doc="Camera instrument", 

68 name="camera", 

69 storageClass="Camera", 

70 dimensions=("instrument",), 

71 isCalibration=True, 

72 ) 

73 

74 skymap = connectionTypes.Input( 

75 doc="Skymap used for tract sharding of output catalog.", 

76 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

77 storageClass="SkyMap", 

78 dimensions=("skymap",), 

79 ) 

80 

81 fgcmLookUpTable = connectionTypes.PrerequisiteInput( 

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

83 "chromatic corrections."), 

84 name="fgcmLookUpTable", 

85 storageClass="Catalog", 

86 dimensions=("instrument",), 

87 deferLoad=True, 

88 ) 

89 

90 fgcmVisitCatalog = connectionTypes.Input( 

91 doc="Catalog of visit information for fgcm", 

92 name="fgcmVisitCatalog", 

93 storageClass="Catalog", 

94 dimensions=("instrument",), 

95 deferLoad=True, 

96 ) 

97 

98 fgcmStandardStars = connectionTypes.Input( 

99 doc="Catalog of standard star data from fgcm fit", 

100 name="fgcm_Cycle{cycleNumber}_StandardStars", 

101 storageClass="SimpleCatalog", 

102 dimensions=("instrument",), 

103 deferLoad=True, 

104 ) 

105 

106 fgcmZeropoints = connectionTypes.Input( 

107 doc="Catalog of zeropoints from fgcm fit", 

108 name="fgcm_Cycle{cycleNumber}_Zeropoints", 

109 storageClass="Catalog", 

110 dimensions=("instrument",), 

111 deferLoad=True, 

112 ) 

113 

114 fgcmAtmosphereParameters = connectionTypes.Input( 

115 doc="Catalog of atmosphere parameters from fgcm fit", 

116 name="fgcm_Cycle{cycleNumber}_AtmosphereParameters", 

117 storageClass="Catalog", 

118 dimensions=("instrument",), 

119 deferLoad=True, 

120 ) 

121 

122 refCat = connectionTypes.PrerequisiteInput( 

123 doc="Reference catalog to use for photometric calibration", 

124 name="cal_ref_cat", 

125 storageClass="SimpleCatalog", 

126 dimensions=("skypix",), 

127 deferLoad=True, 

128 multiple=True, 

129 ) 

130 

131 fgcmPhotoCalib = connectionTypes.Output( 

132 doc=("Per-visit photometric calibrations derived from fgcm calibration. " 

133 "These catalogs use detector id for the id and are sorted for " 

134 "fast lookups of a detector."), 

135 name="fgcmPhotoCalibCatalog", 

136 storageClass="ExposureCatalog", 

137 dimensions=("instrument", "visit",), 

138 multiple=True, 

139 ) 

140 

141 fgcmTransmissionAtmosphere = connectionTypes.Output( 

142 doc="Per-visit atmosphere transmission files produced from fgcm calibration", 

143 name="transmission_atmosphere_fgcm", 

144 storageClass="TransmissionCurve", 

145 dimensions=("instrument", 

146 "visit",), 

147 multiple=True, 

148 ) 

149 

150 fgcmOffsets = connectionTypes.Output( 

151 doc="Per-band offsets computed from doReferenceCalibration", 

152 name="fgcmReferenceCalibrationOffsets", 

153 storageClass="Catalog", 

154 dimensions=("instrument",), 

155 multiple=False, 

156 ) 

157 

158 fgcmTractStars = connectionTypes.Output( 

159 doc="Per-tract fgcm calibrated stars.", 

160 name="fgcm_standard_star", 

161 storageClass="ArrowAstropy", 

162 dimensions=("instrument", "tract", "skymap"), 

163 multiple=True, 

164 ) 

165 

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

167 super().__init__(config=config) 

168 

169 if str(int(config.connections.cycleNumber)) != config.connections.cycleNumber: 

170 raise ValueError("cycleNumber must be of integer format") 

171 

172 if not config.doReferenceCalibration: 

173 self.prerequisiteInputs.remove("refCat") 

174 if not config.doAtmosphereOutput: 

175 self.inputs.remove("fgcmAtmosphereParameters") 

176 if not config.doZeropointOutput: 

177 self.inputs.remove("fgcmZeropoints") 

178 if not config.doReferenceCalibration: 

179 self.outputs.remove("fgcmOffsets") 

180 if not config.doTractStars: 

181 del self.skymap 

182 del self.fgcmTractStars 

183 

184 def getSpatialBoundsConnections(self): 

185 return ("fgcmPhotoCalib",) 

186 

187 

188class FgcmOutputProductsConfig(pipeBase.PipelineTaskConfig, 

189 pipelineConnections=FgcmOutputProductsConnections): 

190 """Config for FgcmOutputProductsTask""" 

191 

192 physicalFilterMap = pexConfig.DictField( 

193 doc="Mapping from 'physicalFilter' to band.", 

194 keytype=str, 

195 itemtype=str, 

196 default={}, 

197 ) 

198 # The following fields refer to calibrating from a reference 

199 # catalog, but in the future this might need to be expanded 

200 doReferenceCalibration = pexConfig.Field( 

201 doc=("Transfer 'absolute' calibration from reference catalog? " 

202 "This afterburner step is unnecessary if reference stars " 

203 "were used in the full fit in FgcmFitCycleTask."), 

204 dtype=bool, 

205 default=False, 

206 ) 

207 doAtmosphereOutput = pexConfig.Field( 

208 doc="Output atmospheres in transmission_atmosphere_fgcm format", 

209 dtype=bool, 

210 default=True, 

211 ) 

212 doZeropointOutput = pexConfig.Field( 

213 doc="Output zeropoints in fgcm_photoCalib format", 

214 dtype=bool, 

215 default=True, 

216 ) 

217 doComposeWcsJacobian = pexConfig.Field( 

218 doc="Compose Jacobian of WCS with fgcm calibration for output photoCalib?", 

219 dtype=bool, 

220 default=True, 

221 ) 

222 doApplyMeanChromaticCorrection = pexConfig.Field( 

223 doc="Apply the mean chromatic correction to the zeropoints?", 

224 dtype=bool, 

225 default=True, 

226 ) 

227 doTractStars = pexConfig.Field( 

228 doc="Output tract-sharded standard stars?", 

229 dtype=bool, 

230 default=True, 

231 # default=False, 

232 ) 

233 photoCal = pexConfig.ConfigurableField( 

234 target=PhotoCalTask, 

235 doc="task to perform 'absolute' calibration", 

236 ) 

237 referencePixelizationNside = pexConfig.Field( 

238 doc="Healpix nside to pixelize catalog to compare to reference catalog", 

239 dtype=int, 

240 default=64, 

241 ) 

242 referencePixelizationMinStars = pexConfig.Field( 

243 doc=("Minimum number of stars per healpix pixel to select for comparison" 

244 "to the specified reference catalog"), 

245 dtype=int, 

246 default=200, 

247 ) 

248 referenceMinMatch = pexConfig.Field( 

249 doc="Minimum number of stars matched to reference catalog to be used in statistics", 

250 dtype=int, 

251 default=50, 

252 ) 

253 referencePixelizationNPixels = pexConfig.Field( 

254 doc=("Number of healpix pixels to sample to do comparison. " 

255 "Doing too many will take a long time and not yield any more " 

256 "precise results because the final number is the median offset " 

257 "(per band) from the set of pixels."), 

258 dtype=int, 

259 default=100, 

260 ) 

261 

262 def setDefaults(self): 

263 pexConfig.Config.setDefaults(self) 

264 

265 # In order to transfer the "absolute" calibration from a reference 

266 # catalog to the relatively calibrated FGCM standard stars (one number 

267 # per band), we use the PhotoCalTask to match stars in a sample of healpix 

268 # pixels. These basic settings ensure that only well-measured, good stars 

269 # from the source and reference catalogs are used for the matching. 

270 

271 # applyColorTerms needs to be False if doReferenceCalibration is False, 

272 # as is the new default after DM-16702 

273 self.photoCal.applyColorTerms = False 

274 self.photoCal.fluxField = 'instFlux' 

275 self.photoCal.magErrFloor = 0.003 

276 self.photoCal.match.referenceSelection.doSignalToNoise = True 

277 self.photoCal.match.referenceSelection.signalToNoise.minimum = 10.0 

278 self.photoCal.match.sourceSelection.doSignalToNoise = True 

279 self.photoCal.match.sourceSelection.signalToNoise.minimum = 10.0 

280 self.photoCal.match.sourceSelection.signalToNoise.fluxField = 'instFlux' 

281 self.photoCal.match.sourceSelection.signalToNoise.errField = 'instFluxErr' 

282 self.photoCal.match.sourceSelection.doFlags = True 

283 self.photoCal.match.sourceSelection.flags.good = [] 

284 self.photoCal.match.sourceSelection.flags.bad = ['flag_badStar'] 

285 self.photoCal.match.sourceSelection.doUnresolved = False 

286 self.photoCal.match.sourceSelection.doRequirePrimary = False 

287 

288 

289class FgcmOutputProductsTask(pipeBase.PipelineTask): 

290 """ 

291 Output products from FGCM global calibration. 

292 """ 

293 

294 ConfigClass = FgcmOutputProductsConfig 

295 _DefaultName = "fgcmOutputProducts" 

296 

297 def __init__(self, **kwargs): 

298 super().__init__(**kwargs) 

299 

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

301 handleDict = {} 

302 handleDict['camera'] = butlerQC.get(inputRefs.camera) 

303 handleDict['fgcmLookUpTable'] = butlerQC.get(inputRefs.fgcmLookUpTable) 

304 handleDict['fgcmVisitCatalog'] = butlerQC.get(inputRefs.fgcmVisitCatalog) 

305 handleDict['fgcmStandardStars'] = butlerQC.get(inputRefs.fgcmStandardStars) 

306 

307 if self.config.doZeropointOutput: 

308 handleDict['fgcmZeropoints'] = butlerQC.get(inputRefs.fgcmZeropoints) 

309 photoCalibRefDict = {photoCalibRef.dataId['visit']: 

310 photoCalibRef for photoCalibRef in outputRefs.fgcmPhotoCalib} 

311 

312 if self.config.doAtmosphereOutput: 

313 handleDict['fgcmAtmosphereParameters'] = butlerQC.get(inputRefs.fgcmAtmosphereParameters) 

314 atmRefDict = {atmRef.dataId['visit']: atmRef for 

315 atmRef in outputRefs.fgcmTransmissionAtmosphere} 

316 

317 if self.config.doReferenceCalibration: 

318 refConfig = LoadReferenceObjectsConfig() 

319 self.refObjLoader = ReferenceObjectLoader(dataIds=[ref.datasetRef.dataId 

320 for ref in inputRefs.refCat], 

321 refCats=butlerQC.get(inputRefs.refCat), 

322 name=self.config.connections.refCat, 

323 log=self.log, 

324 config=refConfig) 

325 else: 

326 self.refObjLoader = None 

327 

328 if self.config.doTractStars: 

329 handleDict['skymap'] = butlerQC.get(inputRefs.skymap) 

330 tractStarRefDict = {tractStarRef.dataId["tract"]: tractStarRef for 

331 tractStarRef in outputRefs.fgcmTractStars} 

332 

333 struct = self.run(handleDict, self.config.physicalFilterMap) 

334 

335 # Output the photoCalib exposure catalogs 

336 if struct.photoCalibCatalogs is not None: 

337 self.log.info("Outputting photoCalib catalogs.") 

338 for visit, expCatalog in struct.photoCalibCatalogs: 

339 butlerQC.put(expCatalog, photoCalibRefDict[visit]) 

340 self.log.info("Done outputting photoCalib catalogs.") 

341 

342 # Output the atmospheres 

343 if struct.atmospheres is not None: 

344 self.log.info("Outputting atmosphere transmission files.") 

345 for visit, atm in struct.atmospheres: 

346 butlerQC.put(atm, atmRefDict[visit]) 

347 self.log.info("Done outputting atmosphere files.") 

348 

349 if self.config.doReferenceCalibration: 

350 # Turn offset into simple catalog for persistence if necessary 

351 schema = afwTable.Schema() 

352 schema.addField('offset', type=np.float64, 

353 doc="Post-process calibration offset (mag)") 

354 offsetCat = afwTable.BaseCatalog(schema) 

355 offsetCat.resize(len(struct.offsets)) 

356 offsetCat['offset'][:] = struct.offsets 

357 

358 butlerQC.put(offsetCat, outputRefs.fgcmOffsets) 

359 

360 if self.config.doTractStars: 

361 self.log.info("Outputting standard stars per-tract.") 

362 for tractId, catalog in struct.tractStars: 

363 butlerQC.put(catalog, tractStarRefDict[tractId]) 

364 

365 return 

366 

367 def run(self, handleDict, physicalFilterMap): 

368 """Run the output products task. 

369 

370 Parameters 

371 ---------- 

372 handleDict : `dict` 

373 All handles are `lsst.daf.butler.DeferredDatasetHandle` 

374 handle dictionary with keys: 

375 

376 ``"camera"`` 

377 Camera object (`lsst.afw.cameraGeom.Camera`) 

378 ``"fgcmLookUpTable"`` 

379 handle for the FGCM look-up table. 

380 ``"fgcmVisitCatalog"`` 

381 handle for visit summary catalog. 

382 ``"fgcmStandardStars"`` 

383 handle for the output standard star catalog. 

384 ``"fgcmZeropoints"`` 

385 handle for the zeropoint data catalog. 

386 ``"fgcmAtmosphereParameters"`` 

387 handle for the atmosphere parameter catalog. 

388 ``"fgcmBuildStarsTableConfig"`` 

389 Config for `lsst.fgcmcal.fgcmBuildStarsTableTask`. 

390 ``"skymap"`` 

391 Skymap for sharding standard stars (optional). 

392 

393 physicalFilterMap : `dict` 

394 Dictionary of mappings from physical filter to FGCM band. 

395 

396 Returns 

397 ------- 

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

399 Output structure with keys: 

400 

401 offsets : `np.ndarray` 

402 Final reference offsets, per band. 

403 atmospheres : `generator` [(`int`, `lsst.afw.image.TransmissionCurve`)] 

404 Generator that returns (visit, transmissionCurve) tuples. 

405 photoCalibCatalogs : `generator` [(`int`, `lsst.afw.table.ExposureCatalog`)] 

406 Generator that returns (visit, exposureCatalog) tuples. 

407 """ 

408 stdCat = handleDict['fgcmStandardStars'].get() 

409 md = stdCat.getMetadata() 

410 bands = md.getArray('BANDS') 

411 

412 if self.config.doReferenceCalibration: 

413 lutCat = handleDict['fgcmLookUpTable'].get() 

414 offsets = self._computeReferenceOffsets(stdCat, lutCat, physicalFilterMap, bands) 

415 else: 

416 offsets = np.zeros(len(bands)) 

417 

418 if self.config.doZeropointOutput: 

419 zptCat = handleDict['fgcmZeropoints'].get() 

420 visitCat = handleDict['fgcmVisitCatalog'].get() 

421 

422 pcgen = self._outputZeropoints(handleDict['camera'], zptCat, visitCat, offsets, bands, 

423 physicalFilterMap) 

424 else: 

425 pcgen = None 

426 

427 if self.config.doAtmosphereOutput: 

428 atmCat = handleDict['fgcmAtmosphereParameters'].get() 

429 atmgen = self._outputAtmospheres(handleDict, atmCat) 

430 else: 

431 atmgen = None 

432 

433 if self.config.doTractStars: 

434 skymap = handleDict['skymap'] 

435 tractstargen = self._outputTractStars(skymap, stdCat) 

436 else: 

437 tractstargen = None 

438 

439 retStruct = pipeBase.Struct(offsets=offsets, 

440 atmospheres=atmgen, 

441 tractStars=tractstargen) 

442 retStruct.photoCalibCatalogs = pcgen 

443 

444 return retStruct 

445 

446 def generateTractOutputProducts(self, handleDict, tract, 

447 visitCat, zptCat, atmCat, stdCat, 

448 fgcmBuildStarsConfig): 

449 """ 

450 Generate the output products for a given tract, as specified in the config. 

451 

452 This method is here to have an alternate entry-point for 

453 FgcmCalibrateTract. 

454 

455 Parameters 

456 ---------- 

457 handleDict : `dict` 

458 All handles are `lsst.daf.butler.DeferredDatasetHandle` 

459 handle dictionary with keys: 

460 

461 ``"camera"`` 

462 Camera object (`lsst.afw.cameraGeom.Camera`) 

463 ``"fgcmLookUpTable"`` 

464 handle for the FGCM look-up table. 

465 tract : `int` 

466 Tract number 

467 visitCat : `lsst.afw.table.BaseCatalog` 

468 FGCM visitCat from `FgcmBuildStarsTask` 

469 zptCat : `lsst.afw.table.BaseCatalog` 

470 FGCM zeropoint catalog from `FgcmFitCycleTask` 

471 atmCat : `lsst.afw.table.BaseCatalog` 

472 FGCM atmosphere parameter catalog from `FgcmFitCycleTask` 

473 stdCat : `lsst.afw.table.SimpleCatalog` 

474 FGCM standard star catalog from `FgcmFitCycleTask` 

475 fgcmBuildStarsConfig : `lsst.fgcmcal.FgcmBuildStarsConfig` 

476 Configuration object from `FgcmBuildStarsTask` 

477 

478 Returns 

479 ------- 

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

481 Output structure with keys: 

482 

483 offsets : `np.ndarray` 

484 Final reference offsets, per band. 

485 atmospheres : `generator` [(`int`, `lsst.afw.image.TransmissionCurve`)] 

486 Generator that returns (visit, transmissionCurve) tuples. 

487 photoCalibCatalogs : `generator` [(`int`, `lsst.afw.table.ExposureCatalog`)] 

488 Generator that returns (visit, exposureCatalog) tuples. 

489 """ 

490 physicalFilterMap = fgcmBuildStarsConfig.physicalFilterMap 

491 

492 md = stdCat.getMetadata() 

493 bands = md.getArray('BANDS') 

494 

495 if self.config.doComposeWcsJacobian and not fgcmBuildStarsConfig.doApplyWcsJacobian: 

496 raise RuntimeError("Cannot compose the WCS jacobian if it hasn't been applied " 

497 "in fgcmBuildStarsTask.") 

498 

499 if not self.config.doComposeWcsJacobian and fgcmBuildStarsConfig.doApplyWcsJacobian: 

500 self.log.warning("Jacobian was applied in build-stars but doComposeWcsJacobian is not set.") 

501 

502 if self.config.doReferenceCalibration: 

503 lutCat = handleDict['fgcmLookUpTable'].get() 

504 offsets = self._computeReferenceOffsets(stdCat, lutCat, bands, physicalFilterMap) 

505 else: 

506 offsets = np.zeros(len(bands)) 

507 

508 if self.config.doZeropointOutput: 

509 pcgen = self._outputZeropoints(handleDict['camera'], zptCat, visitCat, offsets, bands, 

510 physicalFilterMap) 

511 else: 

512 pcgen = None 

513 

514 if self.config.doAtmosphereOutput: 

515 atmgen = self._outputAtmospheres(handleDict, atmCat) 

516 else: 

517 atmgen = None 

518 

519 retStruct = pipeBase.Struct(offsets=offsets, 

520 atmospheres=atmgen) 

521 retStruct.photoCalibCatalogs = pcgen 

522 

523 return retStruct 

524 

525 def _computeReferenceOffsets(self, stdCat, lutCat, physicalFilterMap, bands): 

526 """ 

527 Compute offsets relative to a reference catalog. 

528 

529 This method splits the star catalog into healpix pixels 

530 and computes the calibration transfer for a sample of 

531 these pixels to approximate the 'absolute' calibration 

532 values (on for each band) to apply to transfer the 

533 absolute scale. 

534 

535 Parameters 

536 ---------- 

537 stdCat : `lsst.afw.table.SimpleCatalog` 

538 FGCM standard stars 

539 lutCat : `lsst.afw.table.SimpleCatalog` 

540 FGCM Look-up table 

541 physicalFilterMap : `dict` 

542 Dictionary of mappings from physical filter to FGCM band. 

543 bands : `list` [`str`] 

544 List of band names from FGCM output 

545 Returns 

546 ------- 

547 offsets : `numpy.array` of floats 

548 Per band zeropoint offsets 

549 """ 

550 

551 # Only use stars that are observed in all the bands that were actually used 

552 # This will ensure that we use the same healpix pixels for the absolute 

553 # calibration of each band. 

554 minObs = stdCat['ngood'].min(axis=1) 

555 

556 goodStars = (minObs >= 1) 

557 stdCat = stdCat[goodStars] 

558 

559 self.log.info("Found %d stars with at least 1 good observation in each band" % 

560 (len(stdCat))) 

561 

562 # Associate each band with the appropriate physicalFilter and make 

563 # filterLabels 

564 filterLabels = [] 

565 

566 lutPhysicalFilters = lutCat[0]['physicalFilters'].split(',') 

567 lutStdPhysicalFilters = lutCat[0]['stdPhysicalFilters'].split(',') 

568 physicalFilterMapBands = list(physicalFilterMap.values()) 

569 physicalFilterMapFilters = list(physicalFilterMap.keys()) 

570 for band in bands: 

571 # Find a physical filter associated from the band by doing 

572 # a reverse lookup on the physicalFilterMap dict 

573 physicalFilterMapIndex = physicalFilterMapBands.index(band) 

574 physicalFilter = physicalFilterMapFilters[physicalFilterMapIndex] 

575 # Find the appropriate fgcm standard physicalFilter 

576 lutPhysicalFilterIndex = lutPhysicalFilters.index(physicalFilter) 

577 stdPhysicalFilter = lutStdPhysicalFilters[lutPhysicalFilterIndex] 

578 filterLabels.append(afwImage.FilterLabel(band=band, 

579 physical=stdPhysicalFilter)) 

580 

581 # We have to make a table for each pixel with flux/fluxErr 

582 # This is a temporary table generated for input to the photoCal task. 

583 # These fluxes are not instFlux (they are top-of-the-atmosphere approximate and 

584 # have had chromatic corrections applied to get to the standard system 

585 # specified by the atmosphere/instrumental parameters), nor are they 

586 # in Jansky (since they don't have a proper absolute calibration: the overall 

587 # zeropoint is estimated from the telescope size, etc.) 

588 sourceMapper = afwTable.SchemaMapper(stdCat.schema) 

589 sourceMapper.addMinimalSchema(afwTable.SimpleTable.makeMinimalSchema()) 

590 sourceMapper.editOutputSchema().addField('instFlux', type=np.float64, 

591 doc="instrumental flux (counts)") 

592 sourceMapper.editOutputSchema().addField('instFluxErr', type=np.float64, 

593 doc="instrumental flux error (counts)") 

594 badStarKey = sourceMapper.editOutputSchema().addField('flag_badStar', 

595 type='Flag', 

596 doc="bad flag") 

597 

598 # Split up the stars 

599 # Note that there is an assumption here that the ra/dec coords stored 

600 # on-disk are in radians, and therefore that starObs['coord_ra'] / 

601 # starObs['coord_dec'] return radians when used as an array of numpy float64s. 

602 ipring = hpg.angle_to_pixel( 

603 self.config.referencePixelizationNside, 

604 stdCat['coord_ra'], 

605 stdCat['coord_dec'], 

606 degrees=False, 

607 ) 

608 h, rev = fgcm.fgcmUtilities.histogram_rev_sorted(ipring) 

609 

610 gdpix, = np.where(h >= self.config.referencePixelizationMinStars) 

611 

612 self.log.info("Found %d pixels (nside=%d) with at least %d good stars" % 

613 (gdpix.size, 

614 self.config.referencePixelizationNside, 

615 self.config.referencePixelizationMinStars)) 

616 

617 if gdpix.size < self.config.referencePixelizationNPixels: 

618 self.log.warning("Found fewer good pixels (%d) than preferred in configuration (%d)" % 

619 (gdpix.size, self.config.referencePixelizationNPixels)) 

620 else: 

621 # Sample out the pixels we want to use 

622 gdpix = np.random.choice(gdpix, size=self.config.referencePixelizationNPixels, replace=False) 

623 

624 results = np.zeros(gdpix.size, dtype=[('hpix', 'i4'), 

625 ('nstar', 'i4', len(bands)), 

626 ('nmatch', 'i4', len(bands)), 

627 ('zp', 'f4', len(bands)), 

628 ('zpErr', 'f4', len(bands))]) 

629 results['hpix'] = ipring[rev[rev[gdpix]]] 

630 

631 # We need a boolean index to deal with catalogs... 

632 selected = np.zeros(len(stdCat), dtype=bool) 

633 

634 refFluxFields = [None]*len(bands) 

635 

636 for p_index, pix in enumerate(gdpix): 

637 i1a = rev[rev[pix]: rev[pix + 1]] 

638 

639 # the stdCat afwTable can only be indexed with boolean arrays, 

640 # and not numpy index arrays (see DM-16497). This little trick 

641 # converts the index array into a boolean array 

642 selected[:] = False 

643 selected[i1a] = True 

644 

645 for b_index, filterLabel in enumerate(filterLabels): 

646 struct = self._computeOffsetOneBand(sourceMapper, badStarKey, b_index, 

647 filterLabel, stdCat, 

648 selected, refFluxFields) 

649 results['nstar'][p_index, b_index] = len(i1a) 

650 results['nmatch'][p_index, b_index] = len(struct.arrays.refMag) 

651 results['zp'][p_index, b_index] = struct.zp 

652 results['zpErr'][p_index, b_index] = struct.sigma 

653 

654 # And compute the summary statistics 

655 offsets = np.zeros(len(bands)) 

656 

657 for b_index, band in enumerate(bands): 

658 # make configurable 

659 ok, = np.where(results['nmatch'][:, b_index] >= self.config.referenceMinMatch) 

660 offsets[b_index] = np.median(results['zp'][ok, b_index]) 

661 # use median absolute deviation to estimate Normal sigma 

662 # see https://en.wikipedia.org/wiki/Median_absolute_deviation 

663 madSigma = 1.4826*np.median(np.abs(results['zp'][ok, b_index] - offsets[b_index])) 

664 self.log.info("Reference catalog offset for %s band: %.12f +/- %.12f", 

665 band, offsets[b_index], madSigma) 

666 

667 return offsets 

668 

669 def _computeOffsetOneBand(self, sourceMapper, badStarKey, 

670 b_index, filterLabel, stdCat, selected, refFluxFields): 

671 """ 

672 Compute the zeropoint offset between the fgcm stdCat and the reference 

673 stars for one pixel in one band 

674 

675 Parameters 

676 ---------- 

677 sourceMapper : `lsst.afw.table.SchemaMapper` 

678 Mapper to go from stdCat to calibratable catalog 

679 badStarKey : `lsst.afw.table.Key` 

680 Key for the field with bad stars 

681 b_index : `int` 

682 Index of the band in the star catalog 

683 filterLabel : `lsst.afw.image.FilterLabel` 

684 filterLabel with band and physical filter 

685 stdCat : `lsst.afw.table.SimpleCatalog` 

686 FGCM standard stars 

687 selected : `numpy.array(dtype=bool)` 

688 Boolean array of which stars are in the pixel 

689 refFluxFields : `list` 

690 List of names of flux fields for reference catalog 

691 """ 

692 

693 sourceCat = afwTable.SimpleCatalog(sourceMapper.getOutputSchema()) 

694 sourceCat.reserve(selected.sum()) 

695 sourceCat.extend(stdCat[selected], mapper=sourceMapper) 

696 sourceCat['instFlux'] = 10.**(stdCat['mag_std_noabs'][selected, b_index]/(-2.5)) 

697 sourceCat['instFluxErr'] = (np.log(10.)/2.5)*(stdCat['magErr_std'][selected, b_index] 

698 * sourceCat['instFlux']) 

699 # Make sure we only use stars that have valid measurements 

700 # (This is perhaps redundant with requirements above that the 

701 # stars be observed in all bands, but it can't hurt) 

702 badStar = (stdCat['mag_std_noabs'][selected, b_index] > 90.0) 

703 for rec in sourceCat[badStar]: 

704 rec.set(badStarKey, True) 

705 

706 exposure = afwImage.ExposureF() 

707 exposure.setFilter(filterLabel) 

708 

709 if refFluxFields[b_index] is None: 

710 # Need to find the flux field in the reference catalog 

711 # to work around limitations of DirectMatch in PhotoCal 

712 ctr = stdCat[0].getCoord() 

713 rad = 0.05*lsst.geom.degrees 

714 refDataTest = self.refObjLoader.loadSkyCircle(ctr, rad, filterLabel.bandLabel) 

715 refFluxFields[b_index] = refDataTest.fluxField 

716 

717 # Make a copy of the config so that we can modify it 

718 calConfig = copy.copy(self.config.photoCal.value) 

719 calConfig.match.referenceSelection.signalToNoise.fluxField = refFluxFields[b_index] 

720 calConfig.match.referenceSelection.signalToNoise.errField = refFluxFields[b_index] + 'Err' 

721 calTask = self.config.photoCal.target(refObjLoader=self.refObjLoader, 

722 config=calConfig, 

723 schema=sourceCat.getSchema()) 

724 

725 struct = calTask.run(exposure, sourceCat) 

726 

727 return struct 

728 

729 def _outputZeropoints(self, camera, zptCat, visitCat, offsets, bands, 

730 physicalFilterMap, tract=None): 

731 """Output the zeropoints in fgcm_photoCalib format. 

732 

733 Parameters 

734 ---------- 

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

736 Camera from the butler. 

737 zptCat : `lsst.afw.table.BaseCatalog` 

738 FGCM zeropoint catalog from `FgcmFitCycleTask`. 

739 visitCat : `lsst.afw.table.BaseCatalog` 

740 FGCM visitCat from `FgcmBuildStarsTask`. 

741 offsets : `numpy.array` 

742 Float array of absolute calibration offsets, one for each filter. 

743 bands : `list` [`str`] 

744 List of band names from FGCM output. 

745 physicalFilterMap : `dict` 

746 Dictionary of mappings from physical filter to FGCM band. 

747 tract: `int`, optional 

748 Tract number to output. Default is None (global calibration) 

749 

750 Returns 

751 ------- 

752 photoCalibCatalogs : `generator` [(`int`, `lsst.afw.table.ExposureCatalog`)] 

753 Generator that returns (visit, exposureCatalog) tuples. 

754 """ 

755 # Select visit/ccds where we have a calibration 

756 # This includes ccds where we were able to interpolate from neighboring 

757 # ccds. 

758 cannot_compute = fgcm.fgcmUtilities.zpFlagDict['CANNOT_COMPUTE_ZEROPOINT'] 

759 selected = (((zptCat['fgcmFlag'] & cannot_compute) == 0) 

760 & (zptCat['fgcmZptVar'] > 0.0) 

761 & (zptCat['fgcmZpt'] > FGCM_ILLEGAL_VALUE) 

762 & (zptCat['fgcmfZptCheb'][:, 0] > 0.0)) 

763 

764 # Log warnings for any visit which has no valid zeropoints 

765 badVisits = np.unique(zptCat['visit'][~selected]) 

766 goodVisits = np.unique(zptCat['visit'][selected]) 

767 allBadVisits = badVisits[~np.isin(badVisits, goodVisits)] 

768 for allBadVisit in allBadVisits: 

769 self.log.warning(f'No suitable photoCalib for visit {allBadVisit}') 

770 

771 # Get a mapping from filtername to the offsets 

772 offsetMapping = {} 

773 for f in physicalFilterMap: 

774 # Not every filter in the map will necesarily have a band. 

775 if physicalFilterMap[f] in bands: 

776 offsetMapping[f] = offsets[bands.index(physicalFilterMap[f])] 

777 

778 # Get a mapping from "ccd" to the ccd index used for the scaling 

779 ccdMapping = {} 

780 for ccdIndex, detector in enumerate(camera): 

781 ccdMapping[detector.getId()] = ccdIndex 

782 

783 # And a mapping to get the flat-field scaling values 

784 scalingMapping = {} 

785 for rec in visitCat: 

786 scalingMapping[rec['visit']] = rec['scaling'] 

787 

788 if self.config.doComposeWcsJacobian: 

789 approxPixelAreaFields = computeApproxPixelAreaFields(camera) 

790 

791 # The zptCat is sorted by visit, which is useful 

792 lastVisit = -1 

793 zptVisitCatalog = None 

794 

795 metadata = dafBase.PropertyList() 

796 metadata.add("COMMENT", "Catalog id is detector id, sorted.") 

797 metadata.add("COMMENT", "Only detectors with data have entries.") 

798 

799 for rec in zptCat[selected]: 

800 # Retrieve overall scaling 

801 scaling = scalingMapping[rec['visit']][ccdMapping[rec['detector']]] 

802 

803 # The postCalibrationOffset describe any zeropoint offsets 

804 # to apply after the fgcm calibration. The first part comes 

805 # from the reference catalog match (used in testing). The 

806 # second part comes from the mean chromatic correction 

807 # (if configured). 

808 postCalibrationOffset = offsetMapping[rec['filtername']] 

809 if self.config.doApplyMeanChromaticCorrection: 

810 postCalibrationOffset += rec['fgcmDeltaChrom'] 

811 

812 fgcmSuperStarField = self._getChebyshevBoundedField(rec['fgcmfZptSstarCheb'], 

813 rec['fgcmfZptChebXyMax']) 

814 # Convert from FGCM AB to nJy 

815 fgcmZptField = self._getChebyshevBoundedField((rec['fgcmfZptCheb']*units.AB).to_value(units.nJy), 

816 rec['fgcmfZptChebXyMax'], 

817 offset=postCalibrationOffset, 

818 scaling=scaling) 

819 

820 if self.config.doComposeWcsJacobian: 

821 

822 fgcmField = afwMath.ProductBoundedField([approxPixelAreaFields[rec['detector']], 

823 fgcmSuperStarField, 

824 fgcmZptField]) 

825 else: 

826 # The photoCalib is just the product of the fgcmSuperStarField and the 

827 # fgcmZptField 

828 fgcmField = afwMath.ProductBoundedField([fgcmSuperStarField, fgcmZptField]) 

829 

830 # The "mean" calibration will be set to the center of the ccd for reference 

831 calibCenter = fgcmField.evaluate(fgcmField.getBBox().getCenter()) 

832 calibErr = (np.log(10.0)/2.5)*calibCenter*np.sqrt(rec['fgcmZptVar']) 

833 photoCalib = afwImage.PhotoCalib(calibrationMean=calibCenter, 

834 calibrationErr=calibErr, 

835 calibration=fgcmField, 

836 isConstant=False) 

837 

838 # Return full per-visit exposure catalogs 

839 if rec['visit'] != lastVisit: 

840 # This is a new visit. If the last visit was not -1, yield 

841 # the ExposureCatalog 

842 if lastVisit > -1: 

843 # ensure that the detectors are in sorted order, for fast lookups 

844 zptVisitCatalog.sort() 

845 yield (int(lastVisit), zptVisitCatalog) 

846 else: 

847 # We need to create a new schema 

848 zptExpCatSchema = afwTable.ExposureTable.makeMinimalSchema() 

849 zptExpCatSchema.addField('visit', type='L', doc='Visit number')