LSST Applications g013ef56533+7c9321ec0f,g042eb84c57+c6cfa41bc3,g199a45376c+0ba108daf9,g1fd858c14a+fcad0d0313,g210f2d0738+c0f94c6586,g262e1987ae+a7e710680e,g29ae962dfc+fb55f2edb0,g2ac17093b6+61d6563b1e,g2b1d02342f+df6f932764,g2cef7863aa+aef1011c0b,g2f7ad74990+c0f94c6586,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53cf87ae69,g47891489e3+4316d04fff,g511e8cfd20+baa56acf6c,g53246c7159+8c5ae1fdc5,g54cd7ddccb+fd7ad03fde,g64539dfbff+c0f94c6586,g67b6fd64d1+4316d04fff,g67fd3c3899+c0f94c6586,g6985122a63+4316d04fff,g74acd417e5+ca833bee28,g786e29fd12+668abc6043,g81db2e9a8d+b2ec8e584f,g87389fa792+8856018cbb,g89139ef638+4316d04fff,g8d7436a09f+0a24083b20,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+11eb8fd5b8,gbf99507273+8c5ae1fdc5,gcdda8b9158+e4c84c9d5c,gce8aa8abaa+8c5ae1fdc5,gd7ef33dd92+4316d04fff,gdab6d2f7ff+ca833bee28,ge410e46f29+4316d04fff,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.40
LSST Data Management Base Package
Loading...
Searching...
No Matches
utilities.py
Go to the documentation of this file.
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
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
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
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
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
768 return atmCat
769
770
771def makeStdSchema(nBands):
772 """
773 Make the standard star schema
774
775 Parameters
776 ----------
777 nBands: `int`
778 Number of bands in standard star catalog
779
780 Returns
781 -------
782 stdSchema: `lsst.afw.table.Schema`
783 """
784
785 stdSchema = afwTable.SimpleTable.makeMinimalSchema()
786 stdSchema.addField('isolated_star_id', type='L', doc='ID from isolated star table')
787 stdSchema.addField('ngood', type='ArrayI', doc='Number of good observations',
788 size=nBands)
789 stdSchema.addField('ntotal', type='ArrayI', doc='Number of total observations',
790 size=nBands)
791 stdSchema.addField('mag_std_noabs', type='ArrayF',
792 doc='Standard magnitude (no absolute calibration)',
793 size=nBands)
794 stdSchema.addField('magErr_std', type='ArrayF',
795 doc='Standard magnitude error',
796 size=nBands)
797 stdSchema.addField('npsfcand', type='ArrayI',
798 doc='Number of observations flagged as psf candidates',
799 size=nBands)
800 stdSchema.addField('delta_aper', type='ArrayF',
801 doc='Delta mag (small - large aperture)',
802 size=nBands)
803
804 return stdSchema
805
806
807def makeStdCat(stdSchema, stdStruct, goodBands):
808 """
809 Make the standard star catalog for persistence
810
811 Parameters
812 ----------
813 stdSchema: `lsst.afw.table.Schema`
814 Standard star catalog schema
815 stdStruct: `numpy.ndarray`
816 Standard star structure in FGCM format
817 goodBands: `list`
818 List of good band names used in stdStruct
819
820 Returns
821 -------
822 stdCat: `lsst.afw.table.BaseCatalog`
823 Standard star catalog for persistence
824 """
825
826 stdCat = afwTable.SimpleCatalog(stdSchema)
827 stdCat.resize(stdStruct.size)
828
829 stdCat['id'][:] = stdStruct['FGCM_ID']
830 stdCat['isolated_star_id'][:] = stdStruct['ALTERNATE_ID']
831 stdCat['coord_ra'][:] = stdStruct['RA'] * geom.degrees
832 stdCat['coord_dec'][:] = stdStruct['DEC'] * geom.degrees
833 stdCat['ngood'][:, :] = stdStruct['NGOOD'][:, :]
834 stdCat['ntotal'][:, :] = stdStruct['NTOTAL'][:, :]
835 stdCat['mag_std_noabs'][:, :] = stdStruct['MAG_STD'][:, :]
836 stdCat['magErr_std'][:, :] = stdStruct['MAGERR_STD'][:, :]
837 if 'NPSFCAND' in stdStruct.dtype.names:
838 stdCat['npsfcand'][:, :] = stdStruct['NPSFCAND'][:, :]
839 stdCat['delta_aper'][:, :] = stdStruct['DELTA_APER'][:, :]
840
841 md = PropertyList()
842 md.set("BANDS", list(goodBands))
843 stdCat.setMetadata(md)
844
845 return stdCat
846
847
849 """
850 Compute the radius associated with a CircularApertureFlux or ApFlux field.
851
852 Parameters
853 ----------
854 fluxField : `str`
855 CircularApertureFlux or ApFlux
856
857 Returns
858 -------
859 apertureRadius : `float`
860 Radius of the aperture field, in pixels.
861
862 Raises
863 ------
864 RuntimeError: Raised if flux field is not a CircularApertureFlux,
865 ApFlux, or apFlux.
866 """
867 # TODO: Move this method to more general stack method in DM-25775
868 m = re.search(r'(CircularApertureFlux|ApFlux|apFlux)_(\d+)_(\d+)_', fluxField)
869
870 if m is None:
871 raise RuntimeError(f"Flux field {fluxField} does not correspond to a CircularApertureFlux or ApFlux")
872
873 apertureRadius = float(m.groups()[1]) + float(m.groups()[2])/10.
874
875 return apertureRadius
876
877
878def extractReferenceMags(refStars, bands, filterMap):
879 """
880 Extract reference magnitudes from refStars for given bands and
881 associated filterMap.
882
883 Parameters
884 ----------
885 refStars : `astropy.table.Table` or `lsst.afw.table.BaseCatalog`
886 FGCM reference star catalog.
887 bands : `list`
888 List of bands for calibration.
889 filterMap: `dict`
890 FGCM mapping of filter to band.
891
892 Returns
893 -------
894 refMag : `np.ndarray`
895 nstar x nband array of reference magnitudes.
896 refMagErr : `np.ndarray`
897 nstar x nband array of reference magnitude errors.
898 """
899 hasAstropyMeta = False
900 try:
901 meta = refStars.meta
902 hasAstropyMeta = True
903 except AttributeError:
904 meta = refStars.getMetadata()
905
906 if 'FILTERNAMES' in meta:
907 if hasAstropyMeta:
908 filternames = meta['FILTERNAMES']
909 else:
910 filternames = meta.getArray('FILTERNAMES')
911
912 # The reference catalog that fgcm wants has one entry per band
913 # in the config file
914 refMag = np.zeros((len(refStars), len(bands)),
915 dtype=refStars['refMag'].dtype) + 99.0
916 refMagErr = np.zeros_like(refMag) + 99.0
917 for i, filtername in enumerate(filternames):
918 # We are allowed to run the fit configured so that we do not
919 # use every column in the reference catalog.
920 try:
921 band = filterMap[filtername]
922 except KeyError:
923 continue
924 try:
925 ind = bands.index(band)
926 except ValueError:
927 continue
928
929 refMag[:, ind] = refStars['refMag'][:, i]
930 refMagErr[:, ind] = refStars['refMagErr'][:, i]
931 else:
932 raise RuntimeError("FGCM reference stars missing FILTERNAMES metadata.")
933
934 return refMag, refMagErr
935
936
937def lookupStaticCalibrations(datasetType, registry, quantumDataId, collections):
938 # For static calibrations, we search with a timespan that has unbounded
939 # begin and end; we'll get an error if there's more than one match (because
940 # then it's not static).
941 timespan = Timespan(begin=None, end=None)
942 result = []
943 # First iterate over all of the data IDs for this dataset type that are
944 # consistent with the quantum data ID.
945 for dataId in registry.queryDataIds(datasetType.dimensions, dataId=quantumDataId):
946 # Find the dataset with this data ID using the unbounded timespan.
947 if ref := registry.findDataset(datasetType, dataId, collections=collections, timespan=timespan):
948 result.append(ref)
949 return result
950
951
952def countDetectors(camera, useScienceDetectors):
953 """Count the detectors in the camera.
954
955 This may be limited to just the science detectors.
956
957 Parameters
958 ----------
959 camera : `lsst.afw.cameraGeom.Camera`
960 Camera object.
961 useScienceDetectors : `bool`, optional
962 Limit to just science detectors?
963 """
964 if not useScienceDetectors:
965 return len(camera)
966
967 nDetector = 0
968 for detector in camera:
969 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE:
970 continue
971 nDetector += 1
972
973 return nDetector
Information about a single exposure of an imaging camera.
Definition VisitInfo.h:68
A BoundedField that evaluate the pixel area of a SkyWcs in angular units.
Defines the fields and offsets for a table.
Definition Schema.h:51
Class for storing ordered metadata with comments.
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
translateVisitCatalog(visitCat)
Definition utilities.py:394
computeApertureRadiusFromName(fluxField)
Definition utilities.py:848
makeStdCat(stdSchema, stdStruct, goodBands)
Definition utilities.py:807
makeZptCat(zptSchema, zpStruct)
Definition utilities.py:659
lookupStaticCalibrations(datasetType, registry, quantumDataId, collections)
Definition utilities.py:937
computePixelAreaFieldDetector(detector, visitInfo=None, areaScaling=1.0, approximate=False)
Definition utilities.py:496
makeConfigDict(config, log, camera, maxIter, resetFitParameters, outputZeropoints, lutFilterNames, tract=None, nCore=1, doPlots=False)
Definition utilities.py:50
countDetectors(camera, useScienceDetectors)
Definition utilities.py:952
makeAtmCat(atmSchema, atmStruct)
Definition utilities.py:738
translateFgcmLut(lutCat, physicalFilterMap)
Definition utilities.py:264
computeApproxPixelAreaFields(camera)
Definition utilities.py:537
makeZptSchema(superStarChebyshevSize, zptChebyshevSize)
Definition utilities.py:573
computeReferencePixelScale(camera, useScienceDetectors=False)
Definition utilities.py:471
extractReferenceMags(refStars, bands, filterMap)
Definition utilities.py:878