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 'nExpPerRun': config.nExpPerRun, 

159 'reserveFraction': config.reserveFraction, 

160 'freezeStdAtmosphere': config.freezeStdAtmosphere, 

161 'precomputeSuperStarInitialCycle': config.precomputeSuperStarInitialCycle, 

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

163 'superStarSubCCDChebyshevOrder': config.superStarSubCcdChebyshevOrder, 

164 'superStarSubCCDTriangular': config.superStarSubCcdTriangular, 

165 'superStarSigmaClip': config.superStarSigmaClip, 

166 'superStarPlotCCDResiduals': config.superStarPlotCcdResiduals, 

167 'superStarForceZeroMean': config.superStarForceZeroMean, 

168 'focalPlaneSigmaClip': config.focalPlaneSigmaClip, 

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

170 'ccdGraySubCCDChebyshevOrder': config.ccdGraySubCcdChebyshevOrder, 

171 'ccdGraySubCCDTriangular': config.ccdGraySubCcdTriangular, 

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

173 'ccdGrayFocalPlaneChebyshevOrder': config.ccdGrayFocalPlaneChebyshevOrder, 

174 'ccdGrayFocalPlaneFitMinCcd': config.ccdGrayFocalPlaneFitMinCcd, 

175 'ccdGrayFocalPlaneMaxStars': config.ccdGrayFocalPlaneMaxStars, 

176 'cycleNumber': config.cycleNumber, 

177 'maxIter': maxIter, 

178 'deltaMagBkgOffsetPercentile': config.deltaMagBkgOffsetPercentile, 

179 'deltaMagBkgPerCcd': config.deltaMagBkgPerCcd, 

180 'UTBoundary': config.utBoundary, 

181 'washMJDs': config.washMjds, 

182 'epochMJDs': config.epochMjds, 

183 'coatingMJDs': config.coatingMjds, 

184 'minObsPerBand': config.minObsPerBand, 

185 'latitude': config.latitude, 

186 'defaultCameraOrientation': config.defaultCameraOrientation, 

187 'brightObsGrayMax': config.brightObsGrayMax, 

188 'minStarPerCCD': config.minStarPerCcd, 

189 'minCCDPerExp': config.minCcdPerExp, 

190 'maxCCDGrayErr': config.maxCcdGrayErr, 

191 'minStarPerExp': config.minStarPerExp, 

192 'minExpPerNight': config.minExpPerNight, 

193 'expGrayInitialCut': config.expGrayInitialCut, 

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

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

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

197 'expGrayRecoverCut': config.expGrayRecoverCut, 

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

199 'expGrayErrRecoverCut': config.expGrayErrRecoverCut, 

200 'refStarSnMin': config.refStarSnMin, 

201 'refStarOutlierNSig': config.refStarOutlierNSig, 

202 'applyRefStarColorCuts': config.applyRefStarColorCuts, 

203 'refStarMaxFracUse': config.refStarMaxFracUse, 

204 'useExposureReferenceOffset': config.useExposureReferenceOffset, 

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

206 'starColorCuts': starColorCutList, 

207 'refStarColorCuts': refStarColorCutList, 

208 'aperCorrFitNBins': config.aperCorrFitNBins, 

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

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

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

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

213 'sigFgcmMaxErr': config.sigFgcmMaxErr, 

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

215 'ccdGrayMaxStarErr': config.ccdGrayMaxStarErr, 

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

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

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

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

220 'sigma0Phot': config.sigma0Phot, 

221 'mapLongitudeRef': config.mapLongitudeRef, 

222 'mapNSide': config.mapNSide, 

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

224 'varMinBand': 2, 

225 'useRetrievedPwv': False, 

226 'useNightlyRetrievedPwv': False, 

227 'pwvRetrievalSmoothBlock': 25, 

228 'useQuadraticPwv': config.useQuadraticPwv, 

229 'useRetrievedTauInit': False, 

230 'tauRetrievalMinCCDPerNight': 500, 

231 'modelMagErrors': config.modelMagErrors, 

232 'instrumentParsPerBand': config.instrumentParsPerBand, 

233 'instrumentSlopeMinDeltaT': config.instrumentSlopeMinDeltaT, 

234 'fitMirrorChromaticity': config.fitMirrorChromaticity, 

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

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

237 'autoPhotometricCutNSig': config.autoPhotometricCutNSig, 

238 'autoHighCutNSig': config.autoHighCutNSig, 

239 'deltaAperInnerRadiusArcsec': config.deltaAperInnerRadiusArcsec, 

240 'deltaAperOuterRadiusArcsec': config.deltaAperOuterRadiusArcsec, 

241 'deltaAperFitMinNgoodObs': config.deltaAperFitMinNgoodObs, 

242 'deltaAperFitPerCcdNx': config.deltaAperFitPerCcdNx, 

243 'deltaAperFitPerCcdNy': config.deltaAperFitPerCcdNy, 

244 'deltaAperFitSpatialNside': config.deltaAperFitSpatialNside, 

245 'doComputeDeltaAperExposures': config.doComputeDeltaAperPerVisit, 

246 'doComputeDeltaAperStars': config.doComputeDeltaAperPerStar, 

247 'doComputeDeltaAperMap': config.doComputeDeltaAperMap, 

248 'doComputeDeltaAperPerCcd': config.doComputeDeltaAperPerCcd, 

249 'printOnly': False, 

250 'quietMode': config.quietMode, 

251 'randomSeed': config.randomSeed, 

252 'outputStars': False, 

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

254 'clobber': True, 

255 'useSedLUT': False, 

256 'resetParameters': resetFitParameters, 

257 'doPlots': doPlots, 

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

259 'outputZeropoints': outputZeropoints} 

260 

261 return configDict 

262 

263 

264def translateFgcmLut(lutCat, physicalFilterMap): 

265 """ 

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

267 

268 Parameters 

269 ---------- 

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

271 Catalog describing the FGCM look-up table 

272 physicalFilterMap: `dict` 

273 Physical filter to band mapping 

274 

275 Returns 

276 ------- 

277 fgcmLut: `lsst.fgcm.FgcmLut` 

278 Lookup table for FGCM 

279 lutIndexVals: `numpy.ndarray` 

280 Numpy array with LUT index information for FGCM 

281 lutStd: `numpy.ndarray` 

282 Numpy array with LUT standard throughput values for FGCM 

283 

284 Notes 

285 ----- 

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

287 """ 

288 

289 # first we need the lutIndexVals 

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

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

292 

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

294 # exceptions in the FGCM code. 

295 

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

297 lutFilterNames.size), 

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

299 lutStdFilterNames.size), 

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

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

302 ('PMBELEVATION', 'f8'), 

303 ('LAMBDANORM', 'f8'), 

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

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

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

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

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

309 ('NCCD', 'i4')]) 

310 

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

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

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

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

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

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

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

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

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

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

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

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

323 

324 # now we need the Standard Values 

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

326 ('PWVSTD', 'f8'), 

327 ('O3STD', 'f8'), 

328 ('TAUSTD', 'f8'), 

329 ('ALPHASTD', 'f8'), 

330 ('ZENITHSTD', 'f8'), 

331 ('LAMBDARANGE', 'f8', 2), 

332 ('LAMBDASTEP', 'f8'), 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

359 

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

361 

362 # And the flattened look-up-table 

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

364 ('I1', 'f4')]) 

365 

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

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

368 

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

370 ('D_O3', 'f4'), 

371 ('D_LNTAU', 'f4'), 

372 ('D_ALPHA', 'f4'), 

373 ('D_SECZENITH', 'f4'), 

374 ('D_LNPWV_I1', 'f4'), 

375 ('D_O3_I1', 'f4'), 

376 ('D_LNTAU_I1', 'f4'), 

377 ('D_ALPHA_I1', 'f4'), 

378 ('D_SECZENITH_I1', 'f4')]) 

379 

380 for name in lutDerivFlat.dtype.names: 

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

382 

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

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

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

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

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

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

389 filterToBand=physicalFilterMap) 

390 

391 return fgcmLut, lutIndexVals, lutStd 

392 

393 

394def translateVisitCatalog(visitCat): 

395 """ 

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

397 

398 Parameters 

399 ---------- 

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

401 FGCM visitCat from `lsst.fgcmcal.FgcmBuildStarsTask` 

402 

403 Returns 

404 ------- 

405 fgcmExpInfo: `numpy.ndarray` 

406 Numpy array for visit information for FGCM 

407 

408 Notes 

409 ----- 

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

411 """ 

412 

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

414 

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

416 ('MJD', 'f8'), 

417 ('EXPTIME', 'f8'), 

418 ('PSFFWHM', 'f8'), 

419 ('DELTA_APER', 'f8'), 

420 ('DELTA_APER_DETECTOR', ('f8', nDetector)), 

421 ('SKYBACKGROUND', 'f8'), 

422 ('DEEPFLAG', 'i2'), 

423 ('TELHA', 'f8'), 

424 ('TELRA', 'f8'), 

425 ('TELDEC', 'f8'), 

426 ('TELROT', 'f8'), 

427 ('PMB', 'f8'), 

428 ('FILTERNAME', 'S50')]) 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

445 

446 return fgcmExpInfo 

447 

448 

449def _computeDefaultVisitInfo(): 

450 """ 

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

452 

453 Returns 

454 ------- 

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

456 The visit info object. 

457 """ 

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

459 # since we are looking for relative scales 

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

461 orientation = 0.0*geom.degrees 

462 

463 # Create a temporary visitInfo for input to createInitialSkyWcs 

464 # The orientation does not matter for the area computation 

465 visitInfo = afwImage.VisitInfo(boresightRaDec=boresight, 

466 boresightRotAngle=orientation, 

467 rotType=afwImage.RotType.SKY) 

468 return visitInfo 

469 

470 

471def computeReferencePixelScale(camera, useScienceDetectors=False): 

472 """ 

473 Compute the median pixel scale in the camera 

474 

475 Returns 

476 ------- 

477 pixelScale : `float` 

478 Average pixel scale (arcsecond) over the camera 

479 useScienceDetectors : `bool`, optional 

480 Limit to just science detectors? 

481 """ 

482 visitInfo = _computeDefaultVisitInfo() 

483 

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

485 for i, detector in enumerate(camera): 

486 if useScienceDetectors: 

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

488 continue 

489 wcs = createInitialSkyWcs(visitInfo, detector, False) 

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

491 

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

493 return np.median(pixelScales[ok]) 

494 

495 

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

497 """ 

498 Compute the reference pixel area field for a detector. 

499 

500 Parameters 

501 ---------- 

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

503 The visit info object. 

504 detector : `lsst.afw.cameraGeom.detector` 

505 The detector object. 

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

507 If not supplied, will compute a default visitInfo. 

508 areaScaling : `float`, optional 

509 Area scaling factor. 

510 approximate : `bool`, optional 

511 Compute Chebyshev approximation to pixel area field? 

512 

513 Returns 

514 ------- 

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

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

517 """ 

518 if visitInfo is None: 

519 visitInfo = _computeDefaultVisitInfo() 

520 

521 wcs = createInitialSkyWcs(visitInfo, detector, False) 

522 bbox = detector.getBBox() 

523 

524 pixelAreaField = afwMath.PixelAreaBoundedField( 

525 bbox, 

526 wcs, 

527 unit=geom.arcseconds, 

528 scaling=areaScaling, 

529 ) 

530 

531 if approximate: 

532 pixelAreaField = afwMath.ChebyshevBoundedField.approximate(pixelAreaField) 

533 

534 return pixelAreaField 

535 

536 

537def computeApproxPixelAreaFields(camera): 

538 """ 

539 Compute the approximate pixel area bounded fields from the camera 

540 geometry. 

541 

542 Parameters 

543 ---------- 

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

545 

546 Returns 

547 ------- 

548 approxPixelAreaFields: `dict` 

549 Dictionary of approximate area fields, keyed with detector ID 

550 """ 

551 

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

553 

554 visitInfo = _computeDefaultVisitInfo() 

555 

556 approxPixelAreaFields = {} 

557 

558 for i, detector in enumerate(camera): 

559 key = detector.getId() 

560 

561 approxAreaField = computePixelAreaFieldDetector( 

562 detector, 

563 visitInfo=visitInfo, 

564 areaScaling=areaScaling, 

565 approximate=True, 

566 ) 

567 

568 approxPixelAreaFields[key] = approxAreaField 

569 

570 return approxPixelAreaFields 

571 

572 

573def makeZptSchema(superStarChebyshevSize, zptChebyshevSize): 

574 """ 

575 Make the zeropoint schema 

576 

577 Parameters 

578 ---------- 

579 superStarChebyshevSize: `int` 

580 Length of the superstar chebyshev array 

581 zptChebyshevSize: `int` 

582 Length of the zeropoint chebyshev array 

583 

584 Returns 

585 ------- 

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

587 """ 

588 

589 zptSchema = afwTable.Schema() 

590 

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

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

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

594 '1: Photometric, used in fit; ' 

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

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

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

598 '16: No zeropoint could be determined; ' 

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

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

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

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

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

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

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

606 size=zptChebyshevSize, 

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

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

609 size=superStarChebyshevSize, 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

631 'for 25% bluest stars') 

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

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

634 'for 25% bluest stars') 

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

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

637 'for 25% reddest stars') 

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

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

640 'for 25% reddest stars') 

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

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

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

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

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

646 'at the time of the exposure') 

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

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

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

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

651 '(value set by deltaMagBkgOffsetPercentile) calibration ' 

652 'stars.')) 

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

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

655 

656 return zptSchema 

657 

658 

659def makeZptCat(zptSchema, zpStruct): 

660 """ 

661 Make the zeropoint catalog for persistence 

662 

663 Parameters 

664 ---------- 

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

666 Zeropoint catalog schema 

667 zpStruct: `numpy.ndarray` 

668 Zeropoint structure from fgcm 

669 

670 Returns 

671 ------- 

672 zptCat: `afwTable.BaseCatalog` 

673 Zeropoint catalog for persistence 

674 """ 

675 

676 zptCat = afwTable.BaseCatalog(zptSchema) 

677 zptCat.reserve(zpStruct.size) 

678 

679 for filterName in zpStruct['FILTERNAME']: 

680 rec = zptCat.addNew() 

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

682 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

710 

711 return zptCat 

712 

713 

714def makeAtmSchema(): 

715 """ 

716 Make the atmosphere schema 

717 

718 Returns 

719 ------- 

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

721 """ 

722 

723 atmSchema = afwTable.Schema() 

724 

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

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

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

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

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

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

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

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

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

734 

735 return atmSchema 

736 

737 

738def makeAtmCat(atmSchema, atmStruct): 

739 """ 

740 Make the atmosphere catalog for persistence 

741 

742 Parameters 

743 ---------- 

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

745 Atmosphere catalog schema 

746 atmStruct: `numpy.ndarray` 

747 Atmosphere structure from fgcm 

748 

749 Returns 

750 ------- 

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

752 Atmosphere catalog for persistence 

753 """ 

754 

755 atmCat = afwTable.BaseCatalog(atmSchema) 

756 atmCat.resize(atmStruct.size) 

757 

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

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

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

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

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

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

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

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

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

767