Coverage for python / lsst / fgcmcal / utilities.py: 9%

1# This file is part of fgcmcal. 

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"""Utility functions for fgcmcal. 

22 

23This file contains utility functions that are used by more than one task, 

24and do not need to be part of a task. 

25""" 

26 

27import numpy as np 

28import os 

29import re 

30 

31from lsst.daf.base import PropertyList 

32from lsst.daf.butler import Timespan 

33import lsst.afw.cameraGeom as afwCameraGeom 

34import lsst.afw.table as afwTable 

35import lsst.afw.image as afwImage 

36import lsst.afw.math as afwMath 

37import lsst.geom as geom 

38from lsst.obs.base import createInitialSkyWcs 

39 

40import fgcm 

41 

42 

43FGCM_EXP_FIELD = 'VISIT' 

44FGCM_CCD_FIELD = 'DETECTOR' 

45FGCM_ILLEGAL_VALUE = -9999.0 

46 

47 

48def makeConfigDict(config, log, camera, maxIter, 

49 resetFitParameters, outputZeropoints, 

50 lutFilterNames, tract=None, nCore=1, doPlots=False): 

51 """ 

52 Make the FGCM fit cycle configuration dict 

53 

54 Parameters 

55 ---------- 

56 config : `lsst.fgcmcal.FgcmFitCycleConfig` 

57 Configuration object 

58 log : `logging.Logger` 

59 Log object. 

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

61 Camera from the butler 

62 maxIter : `int` 

63 Maximum number of iterations 

64 resetFitParameters: `bool` 

65 Reset fit parameters before fitting? 

66 outputZeropoints : `bool` 

67 Compute zeropoints for output? 

68 lutFilterNames : array-like, `str` 

69 Array of physical filter names in the LUT. 

70 tract : `int`, optional 

71 Tract number for extending the output file name for debugging. 

72 Default is None. 

73 nCore : `int`, optional 

74 Number of cores to use. 

75 doPlots : `bool`, optional 

76 Make FGCM QA plots? 

77 

78 Returns 

79 ------- 

80 configDict : `dict` 

81 Configuration dictionary for fgcm 

82 """ 

83 # Extract the bands that are _not_ being fit for fgcm configuration 

84 notFitBands = [b for b in config.bands if b not in config.fitBands] 

85 

86 # process the starColorCuts 

87 starColorCutList = [] 

88 for ccut in config.starColorCuts: 

89 if ccut == 'NO_DATA': 

90 # No color cuts to apply. 

91 break 

92 parts = ccut.split(',') 

93 starColorCutList.append([parts[0], parts[1], float(parts[2]), float(parts[3])]) 

94 

95 # process the refStarColorCuts 

96 refStarColorCutList = [] 

97 for ccut in config.refStarColorCuts: 

98 if ccut == 'NO_DATA': 

99 # No color cuts to apply. 

100 break 

101 parts = ccut.split(',') 

102 refStarColorCutList.append([parts[0], parts[1], float(parts[2]), float(parts[3])]) 

103 

104 # TODO: Having direct access to the mirror area from the camera would be 

105 # useful. See DM-16489. 

106 # Mirror area in cm**2 

107 if config.mirrorArea is None: 

108 mirrorArea = np.pi*(camera.telescopeDiameter*100./2.)**2. 

109 else: 

110 # Convert to square cm. 

111 mirrorArea = config.mirrorArea * 100.**2. 

112 

113 if config.cameraGain is None: 

114 # Get approximate average camera gain: 

115 gains = [amp.getGain() for detector in camera for amp in detector.getAmplifiers()] 

116 cameraGain = float(np.median(gains)) 

117 else: 

118 cameraGain = config.cameraGain 

119 

120 # Cut down the filter map to those that are in the LUT 

121 filterToBand = {filterName: config.physicalFilterMap[filterName] for 

122 filterName in lutFilterNames} 

123 

124 if tract is None: 

125 outfileBase = config.outfileBase 

126 else: 

127 outfileBase = '%s-%06d' % (config.outfileBase, tract) 

128 

129 if config.aperCorrPerCcd: 

130 if config.aperCorrUsePsfFwhm: 

131 seeingField = 'PSFFWHM_DETECTOR' 

132 else: 

133 seeingField = 'DELTA_APER_DETECTOR' 

134 seeingSubExposure = True 

135 else: 

136 if config.aperCorrUsePsfFwhm: 

137 seeingField = 'PSFFWHM' 

138 else: 

139 seeingField = 'DELTA_APER' 

140 seeingSubExposure = False 

141 

142 # create a configuration dictionary for fgcmFitCycle 

143 configDict = {'outfileBase': outfileBase, 

144 'logger': log, 

145 'exposureFile': None, 

146 'obsFile': None, 

147 'indexFile': None, 

148 'lutFile': None, 

149 'mirrorArea': mirrorArea, 

150 'cameraGain': cameraGain, 

151 'ccdStartIndex': camera[0].getId(), 

152 'expField': FGCM_EXP_FIELD, 

153 'ccdField': FGCM_CCD_FIELD, 

154 'seeingField': seeingField, 

155 'fwhmField': 'PSFFWHM', 

156 'skyBrightnessField': 'SKYBACKGROUND', 

157 'deepFlag': 'DEEPFLAG', # unused 

158 'bands': list(config.bands), 

159 'fitBands': list(config.fitBands), 

160 'notFitBands': notFitBands, 

161 'requiredBands': list(config.requiredBands), 

162 'filterToBand': filterToBand, 

163 'logLevel': 'INFO', 

164 'nCore': nCore, 

165 'nStarPerRun': config.nStarPerRun, 

166 'nStarPerGrayRun': config.nStarPerGrayRun, 

167 'nObsPerRun': config.nObsPerRun, 

168 'nObsPerGrayRun': config.nObsPerGrayRun, 

169 'nExpPerRun': config.nExpPerRun, 

170 'reserveFraction': config.reserveFraction, 

171 'freezeStdAtmosphere': config.freezeStdAtmosphere, 

172 'precomputeSuperStarInitialCycle': config.precomputeSuperStarInitialCycle, 

173 'superStarSubCCDDict': dict(config.superStarSubCcdDict), 

174 'superStarSubCCDChebyshevOrder': config.superStarSubCcdChebyshevOrder, 

175 'superStarSubCCDTriangular': config.superStarSubCcdTriangular, 

176 'superStarSigmaClip': config.superStarSigmaClip, 

177 'superStarPlotCCDResiduals': config.superStarPlotCcdResiduals, 

178 'superStarForceZeroMean': config.superStarForceZeroMean, 

179 'focalPlaneSigmaClip': config.focalPlaneSigmaClip, 

180 'ccdGraySubCCDDict': dict(config.ccdGraySubCcdDict), 

181 'ccdGraySubCCDChebyshevOrder': config.ccdGraySubCcdChebyshevOrder, 

182 'ccdGraySubCCDTriangular': config.ccdGraySubCcdTriangular, 

183 'ccdGrayFocalPlaneDict': dict(config.ccdGrayFocalPlaneDict), 

184 'ccdGrayFocalPlaneChebyshevOrder': config.ccdGrayFocalPlaneChebyshevOrder, 

185 'ccdGrayFocalPlaneFitMinCcd': config.ccdGrayFocalPlaneFitMinCcd, 

186 'ccdGrayFocalPlaneMaxStars': config.ccdGrayFocalPlaneMaxStars, 

187 'cycleNumber': config.cycleNumber, 

188 'maxIter': maxIter, 

189 'deltaMagBkgOffsetPercentile': config.deltaMagBkgOffsetPercentile, 

190 'deltaMagBkgPerCcd': config.deltaMagBkgPerCcd, 

191 'UTBoundary': config.utBoundary, 

192 'washMJDs': config.washMjds, 

193 'epochMJDs': config.epochMjds, 

194 'coatingMJDs': config.coatingMjds, 

195 'minObsPerBand': config.minObsPerBand, 

196 'latitude': config.latitude, 

197 'defaultCameraOrientation': config.defaultCameraOrientation, 

198 'brightObsGrayMax': config.brightObsGrayMax, 

199 'minStarPerCCD': config.minStarPerCcd, 

200 'minCCDPerExp': config.minCcdPerExp, 

201 'maxCCDGrayErr': config.maxCcdGrayErr, 

202 'minStarPerExp': config.minStarPerExp, 

203 'minExpPerNight': config.minExpPerNight, 

204 'expGrayInitialCut': config.expGrayInitialCut, 

205 'expFwhmCutDict': dict(config.expFwhmCutDict), 

206 'expGrayPhotometricCutDict': dict(config.expGrayPhotometricCutDict), 

207 'expGrayHighCutDict': dict(config.expGrayHighCutDict), 

208 'expGrayRecoverCut': config.expGrayRecoverCut, 

209 'expVarGrayPhotometricCutDict': dict(config.expVarGrayPhotometricCutDict), 

210 'expGrayErrRecoverCut': config.expGrayErrRecoverCut, 

211 'refStarSnMin': config.refStarSnMin, 

212 'refStarOutlierNSig': config.refStarOutlierNSig, 

213 'applyRefStarColorCuts': config.applyRefStarColorCuts, 

214 'refStarMaxFracUse': config.refStarMaxFracUse, 

215 'useExposureReferenceOffset': config.useExposureReferenceOffset, 

216 'illegalValue': FGCM_ILLEGAL_VALUE, # internally used by fgcm. 

217 'starColorCuts': starColorCutList, 

218 'refStarColorCuts': refStarColorCutList, 

219 'aperCorrFitNBins': config.aperCorrFitNBins, 

220 'aperCorrInputSlopeDict': dict(config.aperCorrInputSlopeDict), 

221 'seeingSubExposure': seeingSubExposure, 

222 'sedBoundaryTermDict': config.sedboundaryterms.toDict()['data'], 

223 'sedTermDict': config.sedterms.toDict()['data'], 

224 'colorSplitBands': list(config.colorSplitBands), 

225 'sigFgcmMaxErr': config.sigFgcmMaxErr, 

226 'sigFgcmMaxEGrayDict': dict(config.sigFgcmMaxEGrayDict), 

227 'ccdGrayMaxStarErr': config.ccdGrayMaxStarErr, 

228 'approxThroughputDict': dict(config.approxThroughputDict), 

229 'sigmaCalRange': list(config.sigmaCalRange), 

230 'sigmaCalFitPercentile': list(config.sigmaCalFitPercentile), 

231 'sigmaCalPlotPercentile': list(config.sigmaCalPlotPercentile), 

232 'sigma0Phot': config.sigma0Phot, 

233 'mapLongitudeRef': config.mapLongitudeRef, 

234 'mapNSide': config.mapNSide, 

235 'varNSig': 100.0, # Turn off 'variable star selection' which doesn't work yet 

236 'varMinBand': 2, 

237 'useNightlyRetrievedPwv': False, 

238 'useQuadraticPwv': config.useQuadraticPwv, 

239 'useRetrievedPwv': config.useRetrievedPwv, 

240 'pwvRetrievalSmoothBlock': config.retrievedPwvSmoothingBlock, 

241 'retrievedPwvBands': list(config.retrievedPwvBands), 

242 'useRetrievedTauInit': False, 

243 'tauRetrievalMinCCDPerNight': 500, 

244 'modelMagErrors': config.modelMagErrors, 

245 'instrumentParsPerBand': config.instrumentParsPerBand, 

246 'instrumentSlopeMinDeltaT': config.instrumentSlopeMinDeltaT, 

247 'fitMirrorChromaticity': config.fitMirrorChromaticity, 

248 'fitCCDChromaticityDict': dict(config.fitCcdChromaticityDict), 

249 'useRepeatabilityForExpGrayCutsDict': dict(config.useRepeatabilityForExpGrayCutsDict), 

250 'autoPhotometricCutNSig': config.autoPhotometricCutNSig, 

251 'autoHighCutNSig': config.autoHighCutNSig, 

252 'deltaAperInnerRadiusArcsec': config.deltaAperInnerRadiusArcsec, 

253 'deltaAperOuterRadiusArcsec': config.deltaAperOuterRadiusArcsec, 

254 'deltaAperFitMinNgoodObs': config.deltaAperFitMinNgoodObs, 

255 'deltaAperFitPerCcdNx': config.deltaAperFitPerCcdNx, 

256 'deltaAperFitPerCcdNy': config.deltaAperFitPerCcdNy, 

257 'deltaAperFitSpatialNside': config.deltaAperFitSpatialNside, 

258 'doComputeDeltaAperExposures': config.doComputeDeltaAperPerVisit, 

259 'doComputeDeltaAperStars': config.doComputeDeltaAperPerStar, 

260 'doComputeDeltaAperMap': config.doComputeDeltaAperMap, 

261 'doComputeDeltaAperPerCcd': config.doComputeDeltaAperPerCcd, 

262 'printOnly': False, 

263 'quietMode': config.quietMode, 

264 'randomSeed': config.randomSeed, 

265 'outputStars': False, 

266 'outputPath': os.path.abspath('.'), 

267 'clobber': True, 

268 'useSedLUT': False, 

269 'resetParameters': resetFitParameters, 

270 'doPlots': doPlots, 

271 'outputFgcmcalZpts': True, # when outputting zpts, use fgcmcal format 

272 'outputZeropoints': outputZeropoints} 

273 

274 return configDict 

275 

276 

277def translateFgcmLut(lutCat, physicalFilterMap): 

278 """ 

279 Translate the FGCM look-up-table into an fgcm-compatible object 

280 

281 Parameters 

282 ---------- 

283 lutCat: `lsst.afw.table.BaseCatalog` 

284 Catalog describing the FGCM look-up table 

285 physicalFilterMap: `dict` 

286 Physical filter to band mapping 

287 

288 Returns 

289 ------- 

290 fgcmLut: `lsst.fgcm.FgcmLut` 

291 Lookup table for FGCM 

292 lutIndexVals: `numpy.ndarray` 

293 Numpy array with LUT index information for FGCM 

294 lutStd: `numpy.ndarray` 

295 Numpy array with LUT standard throughput values for FGCM 

296 

297 Notes 

298 ----- 

299 After running this code, it is wise to `del lutCat` to clear the memory. 

300 """ 

301 

302 # first we need the lutIndexVals 

303 lutFilterNames = np.array(lutCat[0]['physicalFilters'].split(','), dtype='U') 

304 lutStdFilterNames = np.array(lutCat[0]['stdPhysicalFilters'].split(','), dtype='U') 

305 

306 # Note that any discrepancies between config values will raise relevant 

307 # exceptions in the FGCM code. 

308 

309 lutIndexVals = np.zeros(1, dtype=[('FILTERNAMES', lutFilterNames.dtype.str, 

310 lutFilterNames.size), 

311 ('STDFILTERNAMES', lutStdFilterNames.dtype.str, 

312 lutStdFilterNames.size), 

313 ('PMB', 'f8', lutCat[0]['pmb'].size), 

314 ('PMBFACTOR', 'f8', lutCat[0]['pmbFactor'].size), 

315 ('PMBELEVATION', 'f8'), 

316 ('LAMBDANORM', 'f8'), 

317 ('PWV', 'f8', lutCat[0]['pwv'].size), 

318 ('O3', 'f8', lutCat[0]['o3'].size), 

319 ('TAU', 'f8', lutCat[0]['tau'].size), 

320 ('ALPHA', 'f8', lutCat[0]['alpha'].size), 

321 ('ZENITH', 'f8', lutCat[0]['zenith'].size), 

322 ('NCCD', 'i4')]) 

323 

324 lutIndexVals['FILTERNAMES'][:] = lutFilterNames 

325 lutIndexVals['STDFILTERNAMES'][:] = lutStdFilterNames 

326 lutIndexVals['PMB'][:] = lutCat[0]['pmb'] 

327 lutIndexVals['PMBFACTOR'][:] = lutCat[0]['pmbFactor'] 

328 lutIndexVals['PMBELEVATION'] = lutCat[0]['pmbElevation'] 

329 lutIndexVals['LAMBDANORM'] = lutCat[0]['lambdaNorm'] 

330 lutIndexVals['PWV'][:] = lutCat[0]['pwv'] 

331 lutIndexVals['O3'][:] = lutCat[0]['o3'] 

332 lutIndexVals['TAU'][:] = lutCat[0]['tau'] 

333 lutIndexVals['ALPHA'][:] = lutCat[0]['alpha'] 

334 lutIndexVals['ZENITH'][:] = lutCat[0]['zenith'] 

335 lutIndexVals['NCCD'] = lutCat[0]['nCcd'] 

336 

337 # now we need the Standard Values 

338 lutStd = np.zeros(1, dtype=[('PMBSTD', 'f8'), 

339 ('PWVSTD', 'f8'), 

340 ('O3STD', 'f8'), 

341 ('TAUSTD', 'f8'), 

342 ('ALPHASTD', 'f8'), 

343 ('ZENITHSTD', 'f8'), 

344 ('LAMBDARANGE', 'f8', 2), 

345 ('LAMBDASTEP', 'f8'), 

346 ('LAMBDASTD', 'f8', lutFilterNames.size), 

347 ('LAMBDASTDFILTER', 'f8', lutStdFilterNames.size), 

348 ('I0STD', 'f8', lutFilterNames.size), 

349 ('I1STD', 'f8', lutFilterNames.size), 

350 ('I10STD', 'f8', lutFilterNames.size), 

351 ('I2STD', 'f8', lutFilterNames.size), 

352 ('LAMBDAB', 'f8', lutFilterNames.size), 

353 ('ATMLAMBDA', 'f8', lutCat[0]['atmLambda'].size), 

354 ('ATMSTDTRANS', 'f8', lutCat[0]['atmStdTrans'].size)]) 

355 lutStd['PMBSTD'] = lutCat[0]['pmbStd'] 

356 lutStd['PWVSTD'] = lutCat[0]['pwvStd'] 

357 lutStd['O3STD'] = lutCat[0]['o3Std'] 

358 lutStd['TAUSTD'] = lutCat[0]['tauStd'] 

359 lutStd['ALPHASTD'] = lutCat[0]['alphaStd'] 

360 lutStd['ZENITHSTD'] = lutCat[0]['zenithStd'] 

361 lutStd['LAMBDARANGE'][:] = lutCat[0]['lambdaRange'][:] 

362 lutStd['LAMBDASTEP'] = lutCat[0]['lambdaStep'] 

363 lutStd['LAMBDASTD'][:] = lutCat[0]['lambdaStd'] 

364 lutStd['LAMBDASTDFILTER'][:] = lutCat[0]['lambdaStdFilter'] 

365 lutStd['I0STD'][:] = lutCat[0]['i0Std'] 

366 lutStd['I1STD'][:] = lutCat[0]['i1Std'] 

367 lutStd['I10STD'][:] = lutCat[0]['i10Std'] 

368 lutStd['I2STD'][:] = lutCat[0]['i2Std'] 

369 lutStd['LAMBDAB'][:] = lutCat[0]['lambdaB'] 

370 lutStd['ATMLAMBDA'][:] = lutCat[0]['atmLambda'][:] 

371 lutStd['ATMSTDTRANS'][:] = lutCat[0]['atmStdTrans'][:] 

372 

373 lutTypes = [row['luttype'] for row in lutCat] 

374 

375 # And the flattened look-up-table 

376 lutFlat = np.zeros(lutCat[0]['lut'].size, dtype=[('I0', 'f4'), 

377 ('I1', 'f4')]) 

378 

379 lutFlat['I0'][:] = lutCat[lutTypes.index('I0')]['lut'][:] 

380 lutFlat['I1'][:] = lutCat[lutTypes.index('I1')]['lut'][:] 

381 

382 lutDerivFlat = np.zeros(lutCat[0]['lut'].size, dtype=[('D_LNPWV', 'f4'), 

383 ('D_O3', 'f4'), 

384 ('D_LNTAU', 'f4'), 

385 ('D_ALPHA', 'f4'), 

386 ('D_SECZENITH', 'f4'), 

387 ('D_LNPWV_I1', 'f4'), 

388 ('D_O3_I1', 'f4'), 

389 ('D_LNTAU_I1', 'f4'), 

390 ('D_ALPHA_I1', 'f4'), 

391 ('D_SECZENITH_I1', 'f4')]) 

392 

393 for name in lutDerivFlat.dtype.names: 

394 lutDerivFlat[name][:] = lutCat[lutTypes.index(name)]['lut'][:] 

395 

396 # The fgcm.FgcmLUT() class copies all the LUT information into special 

397 # shared memory objects that will not blow up the memory usage when used 

398 # with python multiprocessing. Once all the numbers are copied, the 

399 # references to the temporary objects (lutCat, lutFlat, lutDerivFlat) 

400 # will fall out of scope and can be cleaned up by the garbage collector. 

401 fgcmLut = fgcm.FgcmLUT(lutIndexVals, lutFlat, lutDerivFlat, lutStd, 

402 filterToBand=physicalFilterMap) 

403 

404 return fgcmLut, lutIndexVals, lutStd 

405 

406 

407def translateVisitCatalog(visitCat): 

408 """ 

409 Translate the FGCM visit catalog to an fgcm-compatible object 

410 

411 Parameters 

412 ---------- 

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

414 FGCM visitCat from `lsst.fgcmcal.FgcmBuildStarsTask` 

415 

416 Returns 

417 ------- 

418 fgcmExpInfo: `numpy.ndarray` 

419 Numpy array for visit information for FGCM 

420 

421 Notes 

422 ----- 

423 After running this code, it is wise to `del visitCat` to clear the memory. 

424 """ 

425 

426 nDetector = visitCat.schema["deltaAperDetector"].asKey().getSize() 

427 

428 fgcmExpInfo = np.zeros(len(visitCat), dtype=[('VISIT', 'i8'), 

429 ('MJD', 'f8'), 

430 ('EXPTIME', 'f8'), 

431 ('PSFFWHM', 'f8'), 

432 ('PSFFWHM_DETECTOR', ('f8', nDetector)), 

433 ('DELTA_APER', 'f8'), 

434 ('DELTA_APER_DETECTOR', ('f8', nDetector)), 

435 ('SKYBACKGROUND', 'f8'), 

436 ('DEEPFLAG', 'i2'), 

437 ('TELHA', 'f8'), 

438 ('TELRA', 'f8'), 

439 ('TELDEC', 'f8'), 

440 ('TELROT', 'f8'), 

441 ('PMB', 'f8'), 

442 ('FILTERNAME', 'S50')]) 

443 fgcmExpInfo['VISIT'][:] = visitCat['visit'] 

444 fgcmExpInfo['MJD'][:] = visitCat['mjd'] 

445 fgcmExpInfo['EXPTIME'][:] = visitCat['exptime'] 

446 fgcmExpInfo['DEEPFLAG'][:] = visitCat['deepFlag'] 

447 fgcmExpInfo['TELHA'][:] = visitCat['telha'] 

448 fgcmExpInfo['TELRA'][:] = visitCat['telra'] 

449 fgcmExpInfo['TELDEC'][:] = visitCat['teldec'] 

450 fgcmExpInfo['TELROT'][:] = visitCat['telrot'] 

451 fgcmExpInfo['PMB'][:] = visitCat['pmb'] 

452 fgcmExpInfo['PSFFWHM'][:] = visitCat['psfFwhm'] 

453 fgcmExpInfo['PSFFWHM_DETECTOR'][:] = visitCat['psfFwhmDetector'] 

454 fgcmExpInfo['DELTA_APER'][:] = visitCat['deltaAper'] 

455 fgcmExpInfo['DELTA_APER_DETECTOR'][:, :] = visitCat['deltaAperDetector'] 

456 fgcmExpInfo['SKYBACKGROUND'][:] = visitCat['skyBackground'] 

457 # Note that we have to go through asAstropy() to get a string 

458 # array out of an afwTable. This is faster than a row-by-row loop. 

459 fgcmExpInfo['FILTERNAME'][:] = visitCat.asAstropy()['physicalFilter'] 

460 

461 return fgcmExpInfo 

462 

463 

464def _computeDefaultVisitInfo(): 

465 """ 

466 Compute a default visitInfo for use with pixel scale / jacobians. 

467 

468 Returns 

469 ------- 

470 visitInfo : `lsst.afw.image.VisitInfo` 

471 The visit info object. 

472 """ 

473 # Generate fake WCSs centered at 180/0 to avoid the RA=0/360 problem, 

474 # since we are looking for relative scales 

475 boresight = geom.SpherePoint(180.0*geom.degrees, 0.0*geom.degrees) 

476 orientation = 0.0*geom.degrees 

477 

478 # Create a temporary visitInfo for input to createInitialSkyWcs 

479 # The orientation does not matter for the area computation 

480 visitInfo = afwImage.VisitInfo(boresightRaDec=boresight, 

481 boresightRotAngle=orientation, 

482 rotType=afwImage.RotType.SKY) 

483 return visitInfo 

484 

485 

486def computeReferencePixelScale(camera, useScienceDetectors=False): 

487 """ 

488 Compute the median pixel scale in the camera 

489 

490 Returns 

491 ------- 

492 pixelScale : `float` 

493 Average pixel scale (arcsecond) over the camera 

494 useScienceDetectors : `bool`, optional 

495 Limit to just science detectors? 

496 """ 

497 visitInfo = _computeDefaultVisitInfo() 

498 

499 pixelScales = np.zeros(len(camera)) 

500 for i, detector in enumerate(camera): 

501 if useScienceDetectors: 

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

503 continue 

504 wcs = createInitialSkyWcs(visitInfo, detector, False) 

505 pixelScales[i] = wcs.getPixelScale(detector.getBBox().getCenter()).asArcseconds() 

506 

507 ok, = np.where(pixelScales > 0.0) 

508 return np.median(pixelScales[ok]) 

509 

510 

511def computePixelAreaFieldDetector(detector, visitInfo=None, areaScaling=1.0, approximate=False): 

512 """ 

513 Compute the reference pixel area field for a detector. 

514 

515 Parameters 

516 ---------- 

517 visitInfo : `lsst.afw.image.VisitInfo` 

518 The visit info object. 

519 detector : `lsst.afw.cameraGeom.detector` 

520 The detector object. 

521 visitInfo : `lsst.afw.image.VisitInfo`, optional 

522 If not supplied, will compute a default visitInfo. 

523 areaScaling : `float`, optional 

524 Area scaling factor. 

525 approximate : `bool`, optional 

526 Compute Chebyshev approximation to pixel area field? 

527 

528 Returns 

529 ------- 

530 pixelAreaField : `lsst.afw.math.BoundedField` 

531 Bounded field describing the pixel area from the detector model. 

532 """ 

533 if visitInfo is None: 

534 visitInfo = _computeDefaultVisitInfo() 

535 

536 wcs = createInitialSkyWcs(visitInfo, detector, False) 

537 bbox = detector.getBBox() 

538 

539 pixelAreaField = afwMath.PixelAreaBoundedField( 

540 bbox, 

541 wcs, 

542 unit=geom.arcseconds, 

543 scaling=areaScaling, 

544 ) 

545 

546 if approximate: 

547 pixelAreaField = afwMath.ChebyshevBoundedField.approximate(pixelAreaField) 

548 

549 return pixelAreaField 

550 

551 

552def computeApproxPixelAreaFields(camera): 

553 """ 

554 Compute the approximate pixel area bounded fields from the camera 

555 geometry. 

556 

557 Parameters 

558 ---------- 

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

560 

561 Returns 

562 ------- 

563 approxPixelAreaFields: `dict` 

564 Dictionary of approximate area fields, keyed with detector ID 

565 """ 

566 

567 areaScaling = 1. / computeReferencePixelScale(camera)**2. 

568 

569 visitInfo = _computeDefaultVisitInfo() 

570 

571 approxPixelAreaFields = {} 

572 

573 for i, detector in enumerate(camera): 

574 key = detector.getId() 

575 

576 approxAreaField = computePixelAreaFieldDetector( 

577 detector, 

578 visitInfo=visitInfo, 

579 areaScaling=areaScaling, 

580 approximate=True, 

581 ) 

582 

583 approxPixelAreaFields[key] = approxAreaField 

584 

585 return approxPixelAreaFields 

586 

587 

588def makeZptSchema(superStarChebyshevSize, zptChebyshevSize): 

589 """ 

590 Make the zeropoint schema 

591 

592 Parameters 

593 ---------- 

594 superStarChebyshevSize: `int` 

595 Length of the superstar chebyshev array 

596 zptChebyshevSize: `int` 

597 Length of the zeropoint chebyshev array 

598 

599 Returns 

600 ------- 

601 zptSchema: `lsst.afw.table.schema` 

602 """ 

603 

604 zptSchema = afwTable.Schema() 

605 

606 zptSchema.addField('visit', type=np.int64, doc='Visit number') 

607 zptSchema.addField('detector', type=np.int32, doc='Detector ID number') 

608 zptSchema.addField('fgcmFlag', type=np.int32, doc=('FGCM flag value: ' 

609 '1: Photometric, used in fit; ' 

610 '2: Photometric, not used in fit; ' 

611 '4: Non-photometric, on partly photometric night; ' 

612 '8: Non-photometric, on non-photometric night; ' 

613 '16: No zeropoint could be determined; ' 

614 '32: Too few stars for reliable gray computation')) 

615 zptSchema.addField('fgcmZpt', type=np.float64, doc='FGCM zeropoint (center of CCD)') 

616 zptSchema.addField('fgcmZptErr', type=np.float64, 

617 doc='Error on zeropoint, estimated from repeatability + number of obs') 

618 zptSchema.addField('fgcmfZptChebXyMax', type='ArrayD', size=2, 

619 doc='maximum x/maximum y to scale to apply chebyshev parameters') 

620 zptSchema.addField('fgcmfZptCheb', type='ArrayD', 

621 size=zptChebyshevSize, 

622 doc='Chebyshev parameters (flattened) for zeropoint') 

623 zptSchema.addField('fgcmfZptSstarCheb', type='ArrayD', 

624 size=superStarChebyshevSize, 

625 doc='Chebyshev parameters (flattened) for superStarFlat') 

626 zptSchema.addField('fgcmI0', type=np.float64, doc='Integral of the passband') 

627 zptSchema.addField('fgcmI10', type=np.float64, doc='Normalized chromatic integral') 

628 zptSchema.addField('fgcmR0', type=np.float64, 

629 doc='Retrieved i0 integral, estimated from stars (only for flag 1)') 

630 zptSchema.addField('fgcmR10', type=np.float64, 

631 doc='Retrieved i10 integral, estimated from stars (only for flag 1)') 

632 zptSchema.addField('fgcmGry', type=np.float64, 

633 doc='Estimated gray extinction relative to atmospheric solution; ' 

634 'only for fgcmFlag <= 4 (see fgcmFlag) ') 

635 zptSchema.addField('fgcmDeltaChrom', type=np.float64, 

636 doc='Mean chromatic correction for stars in this ccd; ' 

637 'only for fgcmFlag <= 4 (see fgcmFlag)') 

638 zptSchema.addField('fgcmZptVar', type=np.float64, doc='Variance of zeropoint over ccd') 

639 zptSchema.addField('fgcmTilings', type=np.float64, 

640 doc='Number of photometric tilings used for solution for ccd') 

641 zptSchema.addField('fgcmFpGry', type=np.float64, 

642 doc='Average gray extinction over the full focal plane ' 

643 '(same for all ccds in a visit)') 

644 zptSchema.addField('fgcmFpGryBlue', type=np.float64, 

645 doc='Average gray extinction over the full focal plane ' 

646 'for 25% bluest stars') 

647 zptSchema.addField('fgcmFpGryBlueErr', type=np.float64, 

648 doc='Error on Average gray extinction over the full focal plane ' 

649 'for 25% bluest stars') 

650 zptSchema.addField('fgcmFpGryRed', type=np.float64, 

651 doc='Average gray extinction over the full focal plane ' 

652 'for 25% reddest stars') 

653 zptSchema.addField('fgcmFpGryRedErr', type=np.float64, 

654 doc='Error on Average gray extinction over the full focal plane ' 

655 'for 25% reddest stars') 

656 zptSchema.addField('fgcmFpVar', type=np.float64, 

657 doc='Variance of gray extinction over the full focal plane ' 

658 '(same for all ccds in a visit)') 

659 zptSchema.addField('fgcmDust', type=np.float64, 

660 doc='Gray dust extinction from the primary/corrector' 

661 'at the time of the exposure') 

662 zptSchema.addField('fgcmFlat', type=np.float64, doc='Superstarflat illumination correction') 

663 zptSchema.addField('fgcmAperCorr', type=np.float64, doc='Aperture correction estimated by fgcm') 

664 zptSchema.addField('fgcmDeltaMagBkg', type=np.float64, 

665 doc=('Local background correction from brightest percentile ' 

666 '(value set by deltaMagBkgOffsetPercentile) calibration ' 

667 'stars.')) 

668 zptSchema.addField('exptime', type=np.float32, doc='Exposure time') 

669 zptSchema.addField('filtername', type=str, size=30, doc='Filter name') 

670 

671 return zptSchema 

672 

673 

674def makeZptCat(zptSchema, zpStruct): 

675 """ 

676 Make the zeropoint catalog for persistence 

677 

678 Parameters 

679 ---------- 

680 zptSchema: `lsst.afw.table.Schema` 

681 Zeropoint catalog schema 

682 zpStruct: `numpy.ndarray` 

683 Zeropoint structure from fgcm 

684 

685 Returns 

686 ------- 

687 zptCat: `afwTable.BaseCatalog` 

688 Zeropoint catalog for persistence 

689 """ 

690 

691 zptCat = afwTable.BaseCatalog(zptSchema) 

692 zptCat.reserve(zpStruct.size) 

693 

694 for filterName in zpStruct['FILTERNAME']: 

695 rec = zptCat.addNew() 

696 rec['filtername'] = filterName.decode('utf-8') 

697 

698 zptCat['visit'][:] = zpStruct[FGCM_EXP_FIELD] 

699 zptCat['detector'][:] = zpStruct[FGCM_CCD_FIELD] 

700 zptCat['fgcmFlag'][:] = zpStruct['FGCM_FLAG'] 

701 zptCat['fgcmZpt'][:] = zpStruct['FGCM_ZPT'] 

702 zptCat['fgcmZptErr'][:] = zpStruct['FGCM_ZPTERR'] 

703 zptCat['fgcmfZptChebXyMax'][:, :] = zpStruct['FGCM_FZPT_XYMAX'] 

704 zptCat['fgcmfZptCheb'][:, :] = zpStruct['FGCM_FZPT_CHEB'] 

705 zptCat['fgcmfZptSstarCheb'][:, :] = zpStruct['FGCM_FZPT_SSTAR_CHEB'] 

706 zptCat['fgcmI0'][:] = zpStruct['FGCM_I0'] 

707 zptCat['fgcmI10'][:] = zpStruct['FGCM_I10'] 

708 zptCat['fgcmR0'][:] = zpStruct['FGCM_R0'] 

709 zptCat['fgcmR10'][:] = zpStruct['FGCM_R10'] 

710 zptCat['fgcmGry'][:] = zpStruct['FGCM_GRY'] 

711 zptCat['fgcmDeltaChrom'][:] = zpStruct['FGCM_DELTACHROM'] 

712 zptCat['fgcmZptVar'][:] = zpStruct['FGCM_ZPTVAR'] 

713 zptCat['fgcmTilings'][:] = zpStruct['FGCM_TILINGS'] 

714 zptCat['fgcmFpGry'][:] = zpStruct['FGCM_FPGRY'] 

715 zptCat['fgcmFpGryBlue'][:] = zpStruct['FGCM_FPGRY_CSPLIT'][:, 0] 

716 zptCat['fgcmFpGryBlueErr'][:] = zpStruct['FGCM_FPGRY_CSPLITERR'][:, 0] 

717 zptCat['fgcmFpGryRed'][:] = zpStruct['FGCM_FPGRY_CSPLIT'][:, 2] 

718 zptCat['fgcmFpGryRedErr'][:] = zpStruct['FGCM_FPGRY_CSPLITERR'][:, 2] 

719 zptCat['fgcmFpVar'][:] = zpStruct['FGCM_FPVAR'] 

720 zptCat['fgcmDust'][:] = zpStruct['FGCM_DUST'] 

721 zptCat['fgcmFlat'][:] = zpStruct['FGCM_FLAT'] 

722 zptCat['fgcmAperCorr'][:] = zpStruct['FGCM_APERCORR'] 

723 zptCat['fgcmDeltaMagBkg'][:] = zpStruct['FGCM_DELTAMAGBKG'] 

724 zptCat['exptime'][:] = zpStruct['EXPTIME'] 

725 

726 return zptCat 

727 

728 

729def makeAtmSchema(): 

730 """ 

731 Make the atmosphere schema 

732 

733 Returns 

734 ------- 

735 atmSchema: `lsst.afw.table.Schema` 

736 """ 

737 

738 atmSchema = afwTable.Schema() 

739 

740 atmSchema.addField('visit', type=np.int64, doc='Visit number') 

741 atmSchema.addField('pmb', type=np.float64, doc='Barometric pressure (mb)') 

742 atmSchema.addField('pwv', type=np.float64, doc='Water vapor (mm)') 

743 atmSchema.addField('tau', type=np.float64, doc='Aerosol optical depth') 

744 atmSchema.addField('alpha', type=np.float64, doc='Aerosol slope') 

745 atmSchema.addField('o3', type=np.float64, doc='Ozone (dobson)') 

746 atmSchema.addField('secZenith', type=np.float64, doc='Secant(zenith) (~ airmass)') 

747 atmSchema.addField('cTrans', type=np.float64, doc='Transmission correction factor') 

748 atmSchema.addField('lamStd', type=np.float64, doc='Wavelength for transmission correction') 

749 

750 return atmSchema 

751 

752 

753def makeAtmCat(atmSchema, atmStruct): 

754 """ 

755 Make the atmosphere catalog for persistence 

756 

757 Parameters 

758 ---------- 

759 atmSchema: `lsst.afw.table.Schema` 

760 Atmosphere catalog schema 

761 atmStruct: `numpy.ndarray` 

762 Atmosphere structure from fgcm 

763 

764 Returns 

765 ------- 

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

767 Atmosphere catalog for persistence 

768 """ 

769 

770 atmCat = afwTable.BaseCatalog(atmSchema) 

771 atmCat.resize(atmStruct.size) 

772