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

168 'reserveFraction': config.reserveFraction, 

169 'freezeStdAtmosphere': config.freezeStdAtmosphere, 

170 'precomputeSuperStarInitialCycle': config.precomputeSuperStarInitialCycle, 

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

172 'superStarSubCCDChebyshevOrder': config.superStarSubCcdChebyshevOrder, 

173 'superStarSubCCDTriangular': config.superStarSubCcdTriangular, 

174 'superStarSigmaClip': config.superStarSigmaClip, 

175 'superStarPlotCCDResiduals': config.superStarPlotCcdResiduals, 

176 'superStarForceZeroMean': config.superStarForceZeroMean, 

177 'focalPlaneSigmaClip': config.focalPlaneSigmaClip, 

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

179 'ccdGraySubCCDChebyshevOrder': config.ccdGraySubCcdChebyshevOrder, 

180 'ccdGraySubCCDTriangular': config.ccdGraySubCcdTriangular, 

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

182 'ccdGrayFocalPlaneChebyshevOrder': config.ccdGrayFocalPlaneChebyshevOrder, 

183 'ccdGrayFocalPlaneFitMinCcd': config.ccdGrayFocalPlaneFitMinCcd, 

184 'ccdGrayFocalPlaneMaxStars': config.ccdGrayFocalPlaneMaxStars, 

185 'cycleNumber': config.cycleNumber, 

186 'maxIter': maxIter, 

187 'deltaMagBkgOffsetPercentile': config.deltaMagBkgOffsetPercentile, 

188 'deltaMagBkgPerCcd': config.deltaMagBkgPerCcd, 

189 'UTBoundary': config.utBoundary, 

190 'washMJDs': config.washMjds, 

191 'epochMJDs': config.epochMjds, 

192 'coatingMJDs': config.coatingMjds, 

193 'minObsPerBand': config.minObsPerBand, 

194 'latitude': config.latitude, 

195 'defaultCameraOrientation': config.defaultCameraOrientation, 

196 'brightObsGrayMax': config.brightObsGrayMax, 

197 'minStarPerCCD': config.minStarPerCcd, 

198 'minCCDPerExp': config.minCcdPerExp, 

199 'maxCCDGrayErr': config.maxCcdGrayErr, 

200 'minStarPerExp': config.minStarPerExp, 

201 'minExpPerNight': config.minExpPerNight, 

202 'expGrayInitialCut': config.expGrayInitialCut, 

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

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

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

206 'expGrayRecoverCut': config.expGrayRecoverCut, 

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

208 'expGrayErrRecoverCut': config.expGrayErrRecoverCut, 

209 'refStarSnMin': config.refStarSnMin, 

210 'refStarOutlierNSig': config.refStarOutlierNSig, 

211 'applyRefStarColorCuts': config.applyRefStarColorCuts, 

212 'refStarMaxFracUse': config.refStarMaxFracUse, 

213 'useExposureReferenceOffset': config.useExposureReferenceOffset, 

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

215 'starColorCuts': starColorCutList, 

216 'refStarColorCuts': refStarColorCutList, 

217 'aperCorrFitNBins': config.aperCorrFitNBins, 

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

219 'seeingSubExposure': seeingSubExposure, 

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

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

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

223 'sigFgcmMaxErr': config.sigFgcmMaxErr, 

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

225 'ccdGrayMaxStarErr': config.ccdGrayMaxStarErr, 

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

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

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

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

230 'sigma0Phot': config.sigma0Phot, 

231 'mapLongitudeRef': config.mapLongitudeRef, 

232 'mapNSide': config.mapNSide, 

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

234 'varMinBand': 2, 

235 'useRetrievedPwv': False, 

236 'useNightlyRetrievedPwv': False, 

237 'pwvRetrievalSmoothBlock': 25, 

238 'useQuadraticPwv': config.useQuadraticPwv, 

239 'useRetrievedTauInit': False, 

240 'tauRetrievalMinCCDPerNight': 500, 

241 'modelMagErrors': config.modelMagErrors, 

242 'instrumentParsPerBand': config.instrumentParsPerBand, 

243 'instrumentSlopeMinDeltaT': config.instrumentSlopeMinDeltaT, 

244 'fitMirrorChromaticity': config.fitMirrorChromaticity, 

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

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

247 'autoPhotometricCutNSig': config.autoPhotometricCutNSig, 

248 'autoHighCutNSig': config.autoHighCutNSig, 

249 'deltaAperInnerRadiusArcsec': config.deltaAperInnerRadiusArcsec, 

250 'deltaAperOuterRadiusArcsec': config.deltaAperOuterRadiusArcsec, 

251 'deltaAperFitMinNgoodObs': config.deltaAperFitMinNgoodObs, 

252 'deltaAperFitPerCcdNx': config.deltaAperFitPerCcdNx, 

253 'deltaAperFitPerCcdNy': config.deltaAperFitPerCcdNy, 

254 'deltaAperFitSpatialNside': config.deltaAperFitSpatialNside, 

255 'doComputeDeltaAperExposures': config.doComputeDeltaAperPerVisit, 

256 'doComputeDeltaAperStars': config.doComputeDeltaAperPerStar, 

257 'doComputeDeltaAperMap': config.doComputeDeltaAperMap, 

258 'doComputeDeltaAperPerCcd': config.doComputeDeltaAperPerCcd, 

259 'printOnly': False, 

260 'quietMode': config.quietMode, 

261 'randomSeed': config.randomSeed, 

262 'outputStars': False, 

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

264 'clobber': True, 

265 'useSedLUT': False, 

266 'resetParameters': resetFitParameters, 

267 'doPlots': doPlots, 

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

269 'outputZeropoints': outputZeropoints} 

270 

271 return configDict 

272 

273 

274def translateFgcmLut(lutCat, physicalFilterMap): 

275 """ 

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

277 

278 Parameters 

279 ---------- 

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

281 Catalog describing the FGCM look-up table 

282 physicalFilterMap: `dict` 

283 Physical filter to band mapping 

284 

285 Returns 

286 ------- 

287 fgcmLut: `lsst.fgcm.FgcmLut` 

288 Lookup table for FGCM 

289 lutIndexVals: `numpy.ndarray` 

290 Numpy array with LUT index information for FGCM 

291 lutStd: `numpy.ndarray` 

292 Numpy array with LUT standard throughput values for FGCM 

293 

294 Notes 

295 ----- 

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

297 """ 

298 

299 # first we need the lutIndexVals 

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

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

302 

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

304 # exceptions in the FGCM code. 

305 

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

307 lutFilterNames.size), 

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

309 lutStdFilterNames.size), 

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

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

312 ('PMBELEVATION', 'f8'), 

313 ('LAMBDANORM', 'f8'), 

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

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

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

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

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

319 ('NCCD', 'i4')]) 

320 

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

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

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

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

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

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

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

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

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

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

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

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

333 

334 # now we need the Standard Values 

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

336 ('PWVSTD', 'f8'), 

337 ('O3STD', 'f8'), 

338 ('TAUSTD', 'f8'), 

339 ('ALPHASTD', 'f8'), 

340 ('ZENITHSTD', 'f8'), 

341 ('LAMBDARANGE', 'f8', 2), 

342 ('LAMBDASTEP', 'f8'), 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

369 

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

371 

372 # And the flattened look-up-table 

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

374 ('I1', 'f4')]) 

375 

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

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

378 

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

380 ('D_O3', 'f4'), 

381 ('D_LNTAU', 'f4'), 

382 ('D_ALPHA', 'f4'), 

383 ('D_SECZENITH', 'f4'), 

384 ('D_LNPWV_I1', 'f4'), 

385 ('D_O3_I1', 'f4'), 

386 ('D_LNTAU_I1', 'f4'), 

387 ('D_ALPHA_I1', 'f4'), 

388 ('D_SECZENITH_I1', 'f4')]) 

389 

390 for name in lutDerivFlat.dtype.names: 

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

392 

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

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

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

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

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

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

399 filterToBand=physicalFilterMap) 

400 

401 return fgcmLut, lutIndexVals, lutStd 

402 

403 

404def translateVisitCatalog(visitCat): 

405 """ 

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

407 

408 Parameters 

409 ---------- 

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

411 FGCM visitCat from `lsst.fgcmcal.FgcmBuildStarsTask` 

412 

413 Returns 

414 ------- 

415 fgcmExpInfo: `numpy.ndarray` 

416 Numpy array for visit information for FGCM 

417 

418 Notes 

419 ----- 

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

421 """ 

422 

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

424 

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

426 ('MJD', 'f8'), 

427 ('EXPTIME', 'f8'), 

428 ('PSFFWHM', 'f8'), 

429 ('PSFFWHM_DETECTOR', ('f8', nDetector)), 

430 ('DELTA_APER', 'f8'), 

431 ('DELTA_APER_DETECTOR', ('f8', nDetector)), 

432 ('SKYBACKGROUND', 'f8'), 

433 ('DEEPFLAG', 'i2'), 

434 ('TELHA', 'f8'), 

435 ('TELRA', 'f8'), 

436 ('TELDEC', 'f8'), 

437 ('TELROT', 'f8'), 

438 ('PMB', 'f8'), 

439 ('FILTERNAME', 'S50')]) 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

457 

458 return fgcmExpInfo 

459 

460 

461def _computeDefaultVisitInfo(): 

462 """ 

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

464 

465 Returns 

466 ------- 

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

468 The visit info object. 

469 """ 

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

471 # since we are looking for relative scales 

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

473 orientation = 0.0*geom.degrees 

474 

475 # Create a temporary visitInfo for input to createInitialSkyWcs 

476 # The orientation does not matter for the area computation 

477 visitInfo = afwImage.VisitInfo(boresightRaDec=boresight, 

478 boresightRotAngle=orientation, 

479 rotType=afwImage.RotType.SKY) 

480 return visitInfo 

481 

482 

483def computeReferencePixelScale(camera, useScienceDetectors=False): 

484 """ 

485 Compute the median pixel scale in the camera 

486 

487 Returns 

488 ------- 

489 pixelScale : `float` 

490 Average pixel scale (arcsecond) over the camera 

491 useScienceDetectors : `bool`, optional 

492 Limit to just science detectors? 

493 """ 

494 visitInfo = _computeDefaultVisitInfo() 

495 

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

497 for i, detector in enumerate(camera): 

498 if useScienceDetectors: 

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

500 continue 

501 wcs = createInitialSkyWcs(visitInfo, detector, False) 

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

503 

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

505 return np.median(pixelScales[ok]) 

506 

507 

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

509 """ 

510 Compute the reference pixel area field for a detector. 

511 

512 Parameters 

513 ---------- 

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

515 The visit info object. 

516 detector : `lsst.afw.cameraGeom.detector` 

517 The detector object. 

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

519 If not supplied, will compute a default visitInfo. 

520 areaScaling : `float`, optional 

521 Area scaling factor. 

522 approximate : `bool`, optional 

523 Compute Chebyshev approximation to pixel area field? 

524 

525 Returns 

526 ------- 

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

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

529 """ 

530 if visitInfo is None: 

531 visitInfo = _computeDefaultVisitInfo() 

532 

533 wcs = createInitialSkyWcs(visitInfo, detector, False) 

534 bbox = detector.getBBox() 

535 

536 pixelAreaField = afwMath.PixelAreaBoundedField( 

537 bbox, 

538 wcs, 

539 unit=geom.arcseconds, 

540 scaling=areaScaling, 

541 ) 

542 

543 if approximate: 

544 pixelAreaField = afwMath.ChebyshevBoundedField.approximate(pixelAreaField) 

545 

546 return pixelAreaField 

547 

548 

549def computeApproxPixelAreaFields(camera): 

550 """ 

551 Compute the approximate pixel area bounded fields from the camera 

552 geometry. 

553 

554 Parameters 

555 ---------- 

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

557 

558 Returns 

559 ------- 

560 approxPixelAreaFields: `dict` 

561 Dictionary of approximate area fields, keyed with detector ID 

562 """ 

563 

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

565 

566 visitInfo = _computeDefaultVisitInfo() 

567 

568 approxPixelAreaFields = {} 

569 

570 for i, detector in enumerate(camera): 

571 key = detector.getId() 

572 

573 approxAreaField = computePixelAreaFieldDetector( 

574 detector, 

575 visitInfo=visitInfo, 

576 areaScaling=areaScaling, 

577 approximate=True, 

578 ) 

579 

580 approxPixelAreaFields[key] = approxAreaField 

581 

582 return approxPixelAreaFields 

583 

584 

585def makeZptSchema(superStarChebyshevSize, zptChebyshevSize): 

586 """ 

587 Make the zeropoint schema 

588 

589 Parameters 

590 ---------- 

591 superStarChebyshevSize: `int` 

592 Length of the superstar chebyshev array 

593 zptChebyshevSize: `int` 

594 Length of the zeropoint chebyshev array 

595 

596 Returns 

597 ------- 

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

599 """ 

600 

601 zptSchema = afwTable.Schema() 

602 

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

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

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

606 '1: Photometric, used in fit; ' 

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

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

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

610 '16: No zeropoint could be determined; ' 

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

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

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

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

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

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

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

618 size=zptChebyshevSize, 

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

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

621 size=superStarChebyshevSize, 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

643 'for 25% bluest stars') 

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

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

646 'for 25% bluest stars') 

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

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

649 'for 25% reddest stars') 

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

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

652 'for 25% reddest stars') 

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

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

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

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

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

658 'at the time of the exposure') 

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

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

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

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

663 '(value set by deltaMagBkgOffsetPercentile) calibration ' 

664 'stars.')) 

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

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

667 

668 return zptSchema 

669 

670 

671def makeZptCat(zptSchema, zpStruct): 

672 """ 

673 Make the zeropoint catalog for persistence 

674 

675 Parameters 

676 ---------- 

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

678 Zeropoint catalog schema 

679 zpStruct: `numpy.ndarray` 

680 Zeropoint structure from fgcm 

681 

682 Returns 

683 ------- 

684 zptCat: `afwTable.BaseCatalog` 

685 Zeropoint catalog for persistence 

686 """ 

687 

688 zptCat = afwTable.BaseCatalog(zptSchema) 

689 zptCat.reserve(zpStruct.size) 

690 

691 for filterName in zpStruct['FILTERNAME']: 

692 rec = zptCat.addNew() 

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

694 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

722 

723 return zptCat 

724 

725 

726def makeAtmSchema(): 

727 """ 

728 Make the atmosphere schema 

729 

730 Returns 

731 ------- 

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

733 """ 

734 

735 atmSchema = afwTable.Schema() 

736 

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

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

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

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

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

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

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

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

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

746 

747 return atmSchema 

748 

749 

750def makeAtmCat(atmSchema, atmStruct): 

751 """ 

752 Make the atmosphere catalog for persistence 

753 

754 Parameters 

755 ---------- 

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

757 Atmosphere catalog schema 

758 atmStruct: `numpy.ndarray` 

759 Atmosphere structure from fgcm 

760 

761 Returns 

762 ------- 

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

764 Atmosphere catalog for persistence 

765 """ 

766 

767 atmCat = afwTable.BaseCatalog(atmSchema) 

768 atmCat.resize(atmStruct.size) 

769 

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

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

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