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 seeingField = 'DELTA_APER_DETECTOR' 

131 else: 

132 seeingField = 'DELTA_APER' 

133 

134 # create a configuration dictionary for fgcmFitCycle 

135 configDict = {'outfileBase': outfileBase, 

136 'logger': log, 

137 'exposureFile': None, 

138 'obsFile': None, 

139 'indexFile': None, 

140 'lutFile': None, 

141 'mirrorArea': mirrorArea, 

142 'cameraGain': cameraGain, 

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

144 'expField': FGCM_EXP_FIELD, 

145 'ccdField': FGCM_CCD_FIELD, 

146 'seeingField': seeingField, 

147 'fwhmField': 'PSFFWHM', 

148 'skyBrightnessField': 'SKYBACKGROUND', 

149 'deepFlag': 'DEEPFLAG', # unused 

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

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

152 'notFitBands': notFitBands, 

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

154 'filterToBand': filterToBand, 

155 'logLevel': 'INFO', 

156 'nCore': nCore, 

157 'nStarPerRun': config.nStarPerRun, 

158 'nStarPerGrayRun': config.nStarPerGrayRun, 

159 'nExpPerRun': config.nExpPerRun, 

160 'reserveFraction': config.reserveFraction, 

161 'freezeStdAtmosphere': config.freezeStdAtmosphere, 

162 'precomputeSuperStarInitialCycle': config.precomputeSuperStarInitialCycle, 

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

164 'superStarSubCCDChebyshevOrder': config.superStarSubCcdChebyshevOrder, 

165 'superStarSubCCDTriangular': config.superStarSubCcdTriangular, 

166 'superStarSigmaClip': config.superStarSigmaClip, 

167 'superStarPlotCCDResiduals': config.superStarPlotCcdResiduals, 

168 'superStarForceZeroMean': config.superStarForceZeroMean, 

169 'focalPlaneSigmaClip': config.focalPlaneSigmaClip, 

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

171 'ccdGraySubCCDChebyshevOrder': config.ccdGraySubCcdChebyshevOrder, 

172 'ccdGraySubCCDTriangular': config.ccdGraySubCcdTriangular, 

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

174 'ccdGrayFocalPlaneChebyshevOrder': config.ccdGrayFocalPlaneChebyshevOrder, 

175 'ccdGrayFocalPlaneFitMinCcd': config.ccdGrayFocalPlaneFitMinCcd, 

176 'ccdGrayFocalPlaneMaxStars': config.ccdGrayFocalPlaneMaxStars, 

177 'cycleNumber': config.cycleNumber, 

178 'maxIter': maxIter, 

179 'deltaMagBkgOffsetPercentile': config.deltaMagBkgOffsetPercentile, 

180 'deltaMagBkgPerCcd': config.deltaMagBkgPerCcd, 

181 'UTBoundary': config.utBoundary, 

182 'washMJDs': config.washMjds, 

183 'epochMJDs': config.epochMjds, 

184 'coatingMJDs': config.coatingMjds, 

185 'minObsPerBand': config.minObsPerBand, 

186 'latitude': config.latitude, 

187 'defaultCameraOrientation': config.defaultCameraOrientation, 

188 'brightObsGrayMax': config.brightObsGrayMax, 

189 'minStarPerCCD': config.minStarPerCcd, 

190 'minCCDPerExp': config.minCcdPerExp, 

191 'maxCCDGrayErr': config.maxCcdGrayErr, 

192 'minStarPerExp': config.minStarPerExp, 

193 'minExpPerNight': config.minExpPerNight, 

194 'expGrayInitialCut': config.expGrayInitialCut, 

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

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

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

198 'expGrayRecoverCut': config.expGrayRecoverCut, 

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

200 'expGrayErrRecoverCut': config.expGrayErrRecoverCut, 

201 'refStarSnMin': config.refStarSnMin, 

202 'refStarOutlierNSig': config.refStarOutlierNSig, 

203 'applyRefStarColorCuts': config.applyRefStarColorCuts, 

204 'refStarMaxFracUse': config.refStarMaxFracUse, 

205 'useExposureReferenceOffset': config.useExposureReferenceOffset, 

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

207 'starColorCuts': starColorCutList, 

208 'refStarColorCuts': refStarColorCutList, 

209 'aperCorrFitNBins': config.aperCorrFitNBins, 

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

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

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

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

214 'sigFgcmMaxErr': config.sigFgcmMaxErr, 

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

216 'ccdGrayMaxStarErr': config.ccdGrayMaxStarErr, 

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

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

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

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

221 'sigma0Phot': config.sigma0Phot, 

222 'mapLongitudeRef': config.mapLongitudeRef, 

223 'mapNSide': config.mapNSide, 

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

225 'varMinBand': 2, 

226 'useRetrievedPwv': False, 

227 'useNightlyRetrievedPwv': False, 

228 'pwvRetrievalSmoothBlock': 25, 

229 'useQuadraticPwv': config.useQuadraticPwv, 

230 'useRetrievedTauInit': False, 

231 'tauRetrievalMinCCDPerNight': 500, 

232 'modelMagErrors': config.modelMagErrors, 

233 'instrumentParsPerBand': config.instrumentParsPerBand, 

234 'instrumentSlopeMinDeltaT': config.instrumentSlopeMinDeltaT, 

235 'fitMirrorChromaticity': config.fitMirrorChromaticity, 

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

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

238 'autoPhotometricCutNSig': config.autoPhotometricCutNSig, 

239 'autoHighCutNSig': config.autoHighCutNSig, 

240 'deltaAperInnerRadiusArcsec': config.deltaAperInnerRadiusArcsec, 

241 'deltaAperOuterRadiusArcsec': config.deltaAperOuterRadiusArcsec, 

242 'deltaAperFitMinNgoodObs': config.deltaAperFitMinNgoodObs, 

243 'deltaAperFitPerCcdNx': config.deltaAperFitPerCcdNx, 

244 'deltaAperFitPerCcdNy': config.deltaAperFitPerCcdNy, 

245 'deltaAperFitSpatialNside': config.deltaAperFitSpatialNside, 

246 'doComputeDeltaAperExposures': config.doComputeDeltaAperPerVisit, 

247 'doComputeDeltaAperStars': config.doComputeDeltaAperPerStar, 

248 'doComputeDeltaAperMap': config.doComputeDeltaAperMap, 

249 'doComputeDeltaAperPerCcd': config.doComputeDeltaAperPerCcd, 

250 'printOnly': False, 

251 'quietMode': config.quietMode, 

252 'randomSeed': config.randomSeed, 

253 'outputStars': False, 

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

255 'clobber': True, 

256 'useSedLUT': False, 

257 'resetParameters': resetFitParameters, 

258 'doPlots': doPlots, 

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

260 'outputZeropoints': outputZeropoints} 

261 

262 return configDict 

263 

264 

265def translateFgcmLut(lutCat, physicalFilterMap): 

266 """ 

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

268 

269 Parameters 

270 ---------- 

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

272 Catalog describing the FGCM look-up table 

273 physicalFilterMap: `dict` 

274 Physical filter to band mapping 

275 

276 Returns 

277 ------- 

278 fgcmLut: `lsst.fgcm.FgcmLut` 

279 Lookup table for FGCM 

280 lutIndexVals: `numpy.ndarray` 

281 Numpy array with LUT index information for FGCM 

282 lutStd: `numpy.ndarray` 

283 Numpy array with LUT standard throughput values for FGCM 

284 

285 Notes 

286 ----- 

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

288 """ 

289 

290 # first we need the lutIndexVals 

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

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

293 

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

295 # exceptions in the FGCM code. 

296 

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

298 lutFilterNames.size), 

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

300 lutStdFilterNames.size), 

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

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

303 ('PMBELEVATION', 'f8'), 

304 ('LAMBDANORM', 'f8'), 

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

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

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

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

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

310 ('NCCD', 'i4')]) 

311 

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

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

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

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

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

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

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

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

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

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

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

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

324 

325 # now we need the Standard Values 

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

327 ('PWVSTD', 'f8'), 

328 ('O3STD', 'f8'), 

329 ('TAUSTD', 'f8'), 

330 ('ALPHASTD', 'f8'), 

331 ('ZENITHSTD', 'f8'), 

332 ('LAMBDARANGE', 'f8', 2), 

333 ('LAMBDASTEP', 'f8'), 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

360 

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

362 

363 # And the flattened look-up-table 

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

365 ('I1', 'f4')]) 

366 

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

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

369 

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

371 ('D_O3', 'f4'), 

372 ('D_LNTAU', 'f4'), 

373 ('D_ALPHA', 'f4'), 

374 ('D_SECZENITH', 'f4'), 

375 ('D_LNPWV_I1', 'f4'), 

376 ('D_O3_I1', 'f4'), 

377 ('D_LNTAU_I1', 'f4'), 

378 ('D_ALPHA_I1', 'f4'), 

379 ('D_SECZENITH_I1', 'f4')]) 

380 

381 for name in lutDerivFlat.dtype.names: 

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

383 

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

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

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

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

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

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

390 filterToBand=physicalFilterMap) 

391 

392 return fgcmLut, lutIndexVals, lutStd 

393 

394 

395def translateVisitCatalog(visitCat): 

396 """ 

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

398 

399 Parameters 

400 ---------- 

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

402 FGCM visitCat from `lsst.fgcmcal.FgcmBuildStarsTask` 

403 

404 Returns 

405 ------- 

406 fgcmExpInfo: `numpy.ndarray` 

407 Numpy array for visit information for FGCM 

408 

409 Notes 

410 ----- 

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

412 """ 

413 

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

415 

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

417 ('MJD', 'f8'), 

418 ('EXPTIME', 'f8'), 

419 ('PSFFWHM', 'f8'), 

420 ('DELTA_APER', 'f8'), 

421 ('DELTA_APER_DETECTOR', ('f8', nDetector)), 

422 ('SKYBACKGROUND', 'f8'), 

423 ('DEEPFLAG', 'i2'), 

424 ('TELHA', 'f8'), 

425 ('TELRA', 'f8'), 

426 ('TELDEC', 'f8'), 

427 ('TELROT', 'f8'), 

428 ('PMB', 'f8'), 

429 ('FILTERNAME', 'S50')]) 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

446 

447 return fgcmExpInfo 

448 

449 

450def _computeDefaultVisitInfo(): 

451 """ 

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

453 

454 Returns 

455 ------- 

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

457 The visit info object. 

458 """ 

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

460 # since we are looking for relative scales 

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

462 orientation = 0.0*geom.degrees 

463 

464 # Create a temporary visitInfo for input to createInitialSkyWcs 

465 # The orientation does not matter for the area computation 

466 visitInfo = afwImage.VisitInfo(boresightRaDec=boresight, 

467 boresightRotAngle=orientation, 

468 rotType=afwImage.RotType.SKY) 

469 return visitInfo 

470 

471 

472def computeReferencePixelScale(camera, useScienceDetectors=False): 

473 """ 

474 Compute the median pixel scale in the camera 

475 

476 Returns 

477 ------- 

478 pixelScale : `float` 

479 Average pixel scale (arcsecond) over the camera 

480 useScienceDetectors : `bool`, optional 

481 Limit to just science detectors? 

482 """ 

483 visitInfo = _computeDefaultVisitInfo() 

484 

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

486 for i, detector in enumerate(camera): 

487 if useScienceDetectors: 

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

489 continue 

490 wcs = createInitialSkyWcs(visitInfo, detector, False) 

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

492 

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

494 return np.median(pixelScales[ok]) 

495 

496 

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

498 """ 

499 Compute the reference pixel area field for a detector. 

500 

501 Parameters 

502 ---------- 

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

504 The visit info object. 

505 detector : `lsst.afw.cameraGeom.detector` 

506 The detector object. 

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

508 If not supplied, will compute a default visitInfo. 

509 areaScaling : `float`, optional 

510 Area scaling factor. 

511 approximate : `bool`, optional 

512 Compute Chebyshev approximation to pixel area field? 

513 

514 Returns 

515 ------- 

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

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

518 """ 

519 if visitInfo is None: 

520 visitInfo = _computeDefaultVisitInfo() 

521 

522 wcs = createInitialSkyWcs(visitInfo, detector, False) 

523 bbox = detector.getBBox() 

524 

525 pixelAreaField = afwMath.PixelAreaBoundedField( 

526 bbox, 

527 wcs, 

528 unit=geom.arcseconds, 

529 scaling=areaScaling, 

530 ) 

531 

532 if approximate: 

533 pixelAreaField = afwMath.ChebyshevBoundedField.approximate(pixelAreaField) 

534 

535 return pixelAreaField 

536 

537 

538def computeApproxPixelAreaFields(camera): 

539 """ 

540 Compute the approximate pixel area bounded fields from the camera 

541 geometry. 

542 

543 Parameters 

544 ---------- 

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

546 

547 Returns 

548 ------- 

549 approxPixelAreaFields: `dict` 

550 Dictionary of approximate area fields, keyed with detector ID 

551 """ 

552 

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

554 

555 visitInfo = _computeDefaultVisitInfo() 

556 

557 approxPixelAreaFields = {} 

558 

559 for i, detector in enumerate(camera): 

560 key = detector.getId() 

561 

562 approxAreaField = computePixelAreaFieldDetector( 

563 detector, 

564 visitInfo=visitInfo, 

565 areaScaling=areaScaling, 

566 approximate=True, 

567 ) 

568 

569 approxPixelAreaFields[key] = approxAreaField 

570 

571 return approxPixelAreaFields 

572 

573 

574def makeZptSchema(superStarChebyshevSize, zptChebyshevSize): 

575 """ 

576 Make the zeropoint schema 

577 

578 Parameters 

579 ---------- 

580 superStarChebyshevSize: `int` 

581 Length of the superstar chebyshev array 

582 zptChebyshevSize: `int` 

583 Length of the zeropoint chebyshev array 

584 

585 Returns 

586 ------- 

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

588 """ 

589 

590 zptSchema = afwTable.Schema() 

591 

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

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

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

595 '1: Photometric, used in fit; ' 

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

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

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

599 '16: No zeropoint could be determined; ' 

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

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

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

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

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

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

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

607 size=zptChebyshevSize, 

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

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

610 size=superStarChebyshevSize, 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

632 'for 25% bluest stars') 

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

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

635 'for 25% bluest stars') 

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

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

638 'for 25% reddest stars') 

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

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

641 'for 25% reddest stars') 

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

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

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

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

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

647 'at the time of the exposure') 

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

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

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

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

652 '(value set by deltaMagBkgOffsetPercentile) calibration ' 

653 'stars.')) 

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

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

656 

657 return zptSchema 

658 

659 

660def makeZptCat(zptSchema, zpStruct): 

661 """ 

662 Make the zeropoint catalog for persistence 

663 

664 Parameters 

665 ---------- 

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

667 Zeropoint catalog schema 

668 zpStruct: `numpy.ndarray` 

669 Zeropoint structure from fgcm 

670 

671 Returns 

672 ------- 

673 zptCat: `afwTable.BaseCatalog` 

674 Zeropoint catalog for persistence 

675 """ 

676 

677 zptCat = afwTable.BaseCatalog(zptSchema) 

678 zptCat.reserve(zpStruct.size) 

679 

680 for filterName in zpStruct['FILTERNAME']: 

681 rec = zptCat.addNew() 

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

683 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

711 

712 return zptCat 

713 

714 

715def makeAtmSchema(): 

716 """ 

717 Make the atmosphere schema 

718 

719 Returns 

720 ------- 

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

722 """ 

723 

724 atmSchema = afwTable.Schema() 

725 

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

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

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

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

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

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

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

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

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

735 

736 return atmSchema 

737 

738 

739def makeAtmCat(atmSchema, atmStruct): 

740 """ 

741 Make the atmosphere catalog for persistence 

742 

743 Parameters 

744 ---------- 

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

746 Atmosphere catalog schema 

747 atmStruct: `numpy.ndarray` 

748 Atmosphere structure from fgcm 

749 

750 Returns 

751 ------- 

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

753 Atmosphere catalog for persistence 

754 """ 

755 

756 atmCat = afwTable.BaseCatalog(atmSchema) 

757 atmCat.resize(atmStruct.size) 

758 

759 atmCat['visit'][:] = atmStruct['VISIT'] 

760 atmCat['pmb'][:] = atmStruct['PMB'] 

761 atmCat['pwv'][:] = atmStruct['PWV'] 

762 atmCat['tau'][:] = atmStruct['TAU'] 

763 atmCat['alpha'][:] = atmStruct['ALPHA'] 

764 atmCat['o3'][:] = atmStruct['O3'] 

765 atmCat['secZenith'][:] = atmStruct['SECZENITH'] 

766 atmCat['cTrans'][:] = atmStruct['CTRANS'] 

767 atmCat['lamStd'][:] = atmStruct['LAMSTD']