LSST Applications g0f08755f38+82efc23009,g12f32b3c4e+e7bdf1200e,g1653933729+a8ce1bb630,g1a0ca8cf93+50eff2b06f,g28da252d5a+52db39f6a5,g2bbee38e9b+37c5a29d61,g2bc492864f+37c5a29d61,g2cdde0e794+c05ff076ad,g3156d2b45e+41e33cbcdc,g347aa1857d+37c5a29d61,g35bb328faa+a8ce1bb630,g3a166c0a6a+37c5a29d61,g3e281a1b8c+fb992f5633,g414038480c+7f03dfc1b0,g41af890bb2+11b950c980,g5fbc88fb19+17cd334064,g6b1c1869cb+12dd639c9a,g781aacb6e4+a8ce1bb630,g80478fca09+72e9651da0,g82479be7b0+04c31367b4,g858d7b2824+82efc23009,g9125e01d80+a8ce1bb630,g9726552aa6+8047e3811d,ga5288a1d22+e532dc0a0b,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gc28159a63d+37c5a29d61,gcf0d15dbbd+2acd6d4d48,gd7358e8bfb+778a810b6e,gda3e153d99+82efc23009,gda6a2b7d83+2acd6d4d48,gdaeeff99f8+1711a396fd,ge2409df99d+6b12de1076,ge79ae78c31+37c5a29d61,gf0baf85859+d0a5978c5a,gf3967379c6+4954f8c433,gfb92a5be7c+82efc23009,gfec2e1e490+2aaed99252,w.2024.46
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.table as afwTable
34import lsst.afw.image as afwImage
35import lsst.afw.math as afwMath
36import lsst.geom as geom
37from lsst.obs.base import createInitialSkyWcs
38
39import fgcm
40
41
42FGCM_EXP_FIELD = 'VISIT'
43FGCM_CCD_FIELD = 'DETECTOR'
44FGCM_ILLEGAL_VALUE = -9999.0
45
46
47def makeConfigDict(config, log, camera, maxIter,
48 resetFitParameters, outputZeropoints,
49 lutFilterNames, tract=None, nCore=1, doPlots=False):
50 """
51 Make the FGCM fit cycle configuration dict
52
53 Parameters
54 ----------
55 config : `lsst.fgcmcal.FgcmFitCycleConfig`
56 Configuration object
57 log : `logging.Logger`
58 Log object.
59 camera : `lsst.afw.cameraGeom.Camera`
60 Camera from the butler
61 maxIter : `int`
62 Maximum number of iterations
63 resetFitParameters: `bool`
64 Reset fit parameters before fitting?
65 outputZeropoints : `bool`
66 Compute zeropoints for output?
67 lutFilterNames : array-like, `str`
68 Array of physical filter names in the LUT.
69 tract : `int`, optional
70 Tract number for extending the output file name for debugging.
71 Default is None.
72 nCore : `int`, optional
73 Number of cores to use.
74 doPlots : `bool`, optional
75 Make FGCM QA plots?
76
77 Returns
78 -------
79 configDict : `dict`
80 Configuration dictionary for fgcm
81 """
82 # Extract the bands that are _not_ being fit for fgcm configuration
83 notFitBands = [b for b in config.bands if b not in config.fitBands]
84
85 # process the starColorCuts
86 starColorCutList = []
87 for ccut in config.starColorCuts:
88 if ccut == 'NO_DATA':
89 # No color cuts to apply.
90 break
91 parts = ccut.split(',')
92 starColorCutList.append([parts[0], parts[1], float(parts[2]), float(parts[3])])
93
94 # process the refStarColorCuts
95 refStarColorCutList = []
96 for ccut in config.refStarColorCuts:
97 if ccut == 'NO_DATA':
98 # No color cuts to apply.
99 break
100 parts = ccut.split(',')
101 refStarColorCutList.append([parts[0], parts[1], float(parts[2]), float(parts[3])])
102
103 # TODO: Having direct access to the mirror area from the camera would be
104 # useful. See DM-16489.
105 # Mirror area in cm**2
106 if config.mirrorArea is None:
107 mirrorArea = np.pi*(camera.telescopeDiameter*100./2.)**2.
108 else:
109 # Convert to square cm.
110 mirrorArea = config.mirrorArea * 100.**2.
111
112 if config.cameraGain is None:
113 # Get approximate average camera gain:
114 gains = [amp.getGain() for detector in camera for amp in detector.getAmplifiers()]
115 cameraGain = float(np.median(gains))
116 else:
117 cameraGain = config.cameraGain
118
119 # Cut down the filter map to those that are in the LUT
120 filterToBand = {filterName: config.physicalFilterMap[filterName] for
121 filterName in lutFilterNames}
122
123 if tract is None:
124 outfileBase = config.outfileBase
125 else:
126 outfileBase = '%s-%06d' % (config.outfileBase, tract)
127
128 # create a configuration dictionary for fgcmFitCycle
129 configDict = {'outfileBase': outfileBase,
130 'logger': log,
131 'exposureFile': None,
132 'obsFile': None,
133 'indexFile': None,
134 'lutFile': None,
135 'mirrorArea': mirrorArea,
136 'cameraGain': cameraGain,
137 'ccdStartIndex': camera[0].getId(),
138 'expField': FGCM_EXP_FIELD,
139 'ccdField': FGCM_CCD_FIELD,
140 'seeingField': 'DELTA_APER',
141 'fwhmField': 'PSFSIGMA',
142 'skyBrightnessField': 'SKYBACKGROUND',
143 'deepFlag': 'DEEPFLAG', # unused
144 'bands': list(config.bands),
145 'fitBands': list(config.fitBands),
146 'notFitBands': notFitBands,
147 'requiredBands': list(config.requiredBands),
148 'filterToBand': filterToBand,
149 'logLevel': 'INFO',
150 'nCore': nCore,
151 'nStarPerRun': config.nStarPerRun,
152 'nExpPerRun': config.nExpPerRun,
153 'reserveFraction': config.reserveFraction,
154 'freezeStdAtmosphere': config.freezeStdAtmosphere,
155 'precomputeSuperStarInitialCycle': config.precomputeSuperStarInitialCycle,
156 'superStarSubCCDDict': dict(config.superStarSubCcdDict),
157 'superStarSubCCDChebyshevOrder': config.superStarSubCcdChebyshevOrder,
158 'superStarSubCCDTriangular': config.superStarSubCcdTriangular,
159 'superStarSigmaClip': config.superStarSigmaClip,
160 'superStarPlotCCDResiduals': config.superStarPlotCcdResiduals,
161 'focalPlaneSigmaClip': config.focalPlaneSigmaClip,
162 'ccdGraySubCCDDict': dict(config.ccdGraySubCcdDict),
163 'ccdGraySubCCDChebyshevOrder': config.ccdGraySubCcdChebyshevOrder,
164 'ccdGraySubCCDTriangular': config.ccdGraySubCcdTriangular,
165 'ccdGrayFocalPlaneDict': dict(config.ccdGrayFocalPlaneDict),
166 'ccdGrayFocalPlaneChebyshevOrder': config.ccdGrayFocalPlaneChebyshevOrder,
167 'ccdGrayFocalPlaneFitMinCcd': config.ccdGrayFocalPlaneFitMinCcd,
168 'cycleNumber': config.cycleNumber,
169 'maxIter': maxIter,
170 'deltaMagBkgOffsetPercentile': config.deltaMagBkgOffsetPercentile,
171 'deltaMagBkgPerCcd': config.deltaMagBkgPerCcd,
172 'UTBoundary': config.utBoundary,
173 'washMJDs': config.washMjds,
174 'epochMJDs': config.epochMjds,
175 'coatingMJDs': config.coatingMjds,
176 'minObsPerBand': config.minObsPerBand,
177 'latitude': config.latitude,
178 'defaultCameraOrientation': config.defaultCameraOrientation,
179 'brightObsGrayMax': config.brightObsGrayMax,
180 'minStarPerCCD': config.minStarPerCcd,
181 'minCCDPerExp': config.minCcdPerExp,
182 'maxCCDGrayErr': config.maxCcdGrayErr,
183 'minStarPerExp': config.minStarPerExp,
184 'minExpPerNight': config.minExpPerNight,
185 'expGrayInitialCut': config.expGrayInitialCut,
186 'expGrayPhotometricCutDict': dict(config.expGrayPhotometricCutDict),
187 'expGrayHighCutDict': dict(config.expGrayHighCutDict),
188 'expGrayRecoverCut': config.expGrayRecoverCut,
189 'expVarGrayPhotometricCutDict': dict(config.expVarGrayPhotometricCutDict),
190 'expGrayErrRecoverCut': config.expGrayErrRecoverCut,
191 'refStarSnMin': config.refStarSnMin,
192 'refStarOutlierNSig': config.refStarOutlierNSig,
193 'applyRefStarColorCuts': config.applyRefStarColorCuts,
194 'refStarMaxFracUse': config.refStarMaxFracUse,
195 'useExposureReferenceOffset': config.useExposureReferenceOffset,
196 'illegalValue': FGCM_ILLEGAL_VALUE, # internally used by fgcm.
197 'starColorCuts': starColorCutList,
198 'refStarColorCuts': refStarColorCutList,
199 'aperCorrFitNBins': config.aperCorrFitNBins,
200 'aperCorrInputSlopeDict': dict(config.aperCorrInputSlopeDict),
201 'sedBoundaryTermDict': config.sedboundaryterms.toDict()['data'],
202 'sedTermDict': config.sedterms.toDict()['data'],
203 'colorSplitBands': list(config.colorSplitBands),
204 'sigFgcmMaxErr': config.sigFgcmMaxErr,
205 'sigFgcmMaxEGrayDict': dict(config.sigFgcmMaxEGrayDict),
206 'ccdGrayMaxStarErr': config.ccdGrayMaxStarErr,
207 'approxThroughputDict': dict(config.approxThroughputDict),
208 'sigmaCalRange': list(config.sigmaCalRange),
209 'sigmaCalFitPercentile': list(config.sigmaCalFitPercentile),
210 'sigmaCalPlotPercentile': list(config.sigmaCalPlotPercentile),
211 'sigma0Phot': config.sigma0Phot,
212 'mapLongitudeRef': config.mapLongitudeRef,
213 'mapNSide': config.mapNSide,
214 'varNSig': 100.0, # Turn off 'variable star selection' which doesn't work yet
215 'varMinBand': 2,
216 'useRetrievedPwv': False,
217 'useNightlyRetrievedPwv': False,
218 'pwvRetrievalSmoothBlock': 25,
219 'useQuadraticPwv': config.useQuadraticPwv,
220 'useRetrievedTauInit': False,
221 'tauRetrievalMinCCDPerNight': 500,
222 'modelMagErrors': config.modelMagErrors,
223 'instrumentParsPerBand': config.instrumentParsPerBand,
224 'instrumentSlopeMinDeltaT': config.instrumentSlopeMinDeltaT,
225 'fitMirrorChromaticity': config.fitMirrorChromaticity,
226 'fitCCDChromaticityDict': dict(config.fitCcdChromaticityDict),
227 'useRepeatabilityForExpGrayCutsDict': dict(config.useRepeatabilityForExpGrayCutsDict),
228 'autoPhotometricCutNSig': config.autoPhotometricCutNSig,
229 'autoHighCutNSig': config.autoHighCutNSig,
230 'deltaAperInnerRadiusArcsec': config.deltaAperInnerRadiusArcsec,
231 'deltaAperOuterRadiusArcsec': config.deltaAperOuterRadiusArcsec,
232 'deltaAperFitMinNgoodObs': config.deltaAperFitMinNgoodObs,
233 'deltaAperFitPerCcdNx': config.deltaAperFitPerCcdNx,
234 'deltaAperFitPerCcdNy': config.deltaAperFitPerCcdNy,
235 'deltaAperFitSpatialNside': config.deltaAperFitSpatialNside,
236 'doComputeDeltaAperExposures': config.doComputeDeltaAperPerVisit,
237 'doComputeDeltaAperStars': config.doComputeDeltaAperPerStar,
238 'doComputeDeltaAperMap': config.doComputeDeltaAperMap,
239 'doComputeDeltaAperPerCcd': config.doComputeDeltaAperPerCcd,
240 'printOnly': False,
241 'quietMode': config.quietMode,
242 'randomSeed': config.randomSeed,
243 'outputStars': False,
244 'outputPath': os.path.abspath('.'),
245 'clobber': True,
246 'useSedLUT': False,
247 'resetParameters': resetFitParameters,
248 'doPlots': doPlots,
249 'outputFgcmcalZpts': True, # when outputting zpts, use fgcmcal format
250 'outputZeropoints': outputZeropoints}
251
252 return configDict
253
254
255def translateFgcmLut(lutCat, physicalFilterMap):
256 """
257 Translate the FGCM look-up-table into an fgcm-compatible object
258
259 Parameters
260 ----------
261 lutCat: `lsst.afw.table.BaseCatalog`
262 Catalog describing the FGCM look-up table
263 physicalFilterMap: `dict`
264 Physical filter to band mapping
265
266 Returns
267 -------
268 fgcmLut: `lsst.fgcm.FgcmLut`
269 Lookup table for FGCM
270 lutIndexVals: `numpy.ndarray`
271 Numpy array with LUT index information for FGCM
272 lutStd: `numpy.ndarray`
273 Numpy array with LUT standard throughput values for FGCM
274
275 Notes
276 -----
277 After running this code, it is wise to `del lutCat` to clear the memory.
278 """
279
280 # first we need the lutIndexVals
281 lutFilterNames = np.array(lutCat[0]['physicalFilters'].split(','), dtype='U')
282 lutStdFilterNames = np.array(lutCat[0]['stdPhysicalFilters'].split(','), dtype='U')
283
284 # Note that any discrepancies between config values will raise relevant
285 # exceptions in the FGCM code.
286
287 lutIndexVals = np.zeros(1, dtype=[('FILTERNAMES', lutFilterNames.dtype.str,
288 lutFilterNames.size),
289 ('STDFILTERNAMES', lutStdFilterNames.dtype.str,
290 lutStdFilterNames.size),
291 ('PMB', 'f8', lutCat[0]['pmb'].size),
292 ('PMBFACTOR', 'f8', lutCat[0]['pmbFactor'].size),
293 ('PMBELEVATION', 'f8'),
294 ('LAMBDANORM', 'f8'),
295 ('PWV', 'f8', lutCat[0]['pwv'].size),
296 ('O3', 'f8', lutCat[0]['o3'].size),
297 ('TAU', 'f8', lutCat[0]['tau'].size),
298 ('ALPHA', 'f8', lutCat[0]['alpha'].size),
299 ('ZENITH', 'f8', lutCat[0]['zenith'].size),
300 ('NCCD', 'i4')])
301
302 lutIndexVals['FILTERNAMES'][:] = lutFilterNames
303 lutIndexVals['STDFILTERNAMES'][:] = lutStdFilterNames
304 lutIndexVals['PMB'][:] = lutCat[0]['pmb']
305 lutIndexVals['PMBFACTOR'][:] = lutCat[0]['pmbFactor']
306 lutIndexVals['PMBELEVATION'] = lutCat[0]['pmbElevation']
307 lutIndexVals['LAMBDANORM'] = lutCat[0]['lambdaNorm']
308 lutIndexVals['PWV'][:] = lutCat[0]['pwv']
309 lutIndexVals['O3'][:] = lutCat[0]['o3']
310 lutIndexVals['TAU'][:] = lutCat[0]['tau']
311 lutIndexVals['ALPHA'][:] = lutCat[0]['alpha']
312 lutIndexVals['ZENITH'][:] = lutCat[0]['zenith']
313 lutIndexVals['NCCD'] = lutCat[0]['nCcd']
314
315 # now we need the Standard Values
316 lutStd = np.zeros(1, dtype=[('PMBSTD', 'f8'),
317 ('PWVSTD', 'f8'),
318 ('O3STD', 'f8'),
319 ('TAUSTD', 'f8'),
320 ('ALPHASTD', 'f8'),
321 ('ZENITHSTD', 'f8'),
322 ('LAMBDARANGE', 'f8', 2),
323 ('LAMBDASTEP', 'f8'),
324 ('LAMBDASTD', 'f8', lutFilterNames.size),
325 ('LAMBDASTDFILTER', 'f8', lutStdFilterNames.size),
326 ('I0STD', 'f8', lutFilterNames.size),
327 ('I1STD', 'f8', lutFilterNames.size),
328 ('I10STD', 'f8', lutFilterNames.size),
329 ('I2STD', 'f8', lutFilterNames.size),
330 ('LAMBDAB', 'f8', lutFilterNames.size),
331 ('ATMLAMBDA', 'f8', lutCat[0]['atmLambda'].size),
332 ('ATMSTDTRANS', 'f8', lutCat[0]['atmStdTrans'].size)])
333 lutStd['PMBSTD'] = lutCat[0]['pmbStd']
334 lutStd['PWVSTD'] = lutCat[0]['pwvStd']
335 lutStd['O3STD'] = lutCat[0]['o3Std']
336 lutStd['TAUSTD'] = lutCat[0]['tauStd']
337 lutStd['ALPHASTD'] = lutCat[0]['alphaStd']
338 lutStd['ZENITHSTD'] = lutCat[0]['zenithStd']
339 lutStd['LAMBDARANGE'][:] = lutCat[0]['lambdaRange'][:]
340 lutStd['LAMBDASTEP'] = lutCat[0]['lambdaStep']
341 lutStd['LAMBDASTD'][:] = lutCat[0]['lambdaStd']
342 lutStd['LAMBDASTDFILTER'][:] = lutCat[0]['lambdaStdFilter']
343 lutStd['I0STD'][:] = lutCat[0]['i0Std']
344 lutStd['I1STD'][:] = lutCat[0]['i1Std']
345 lutStd['I10STD'][:] = lutCat[0]['i10Std']
346 lutStd['I2STD'][:] = lutCat[0]['i2Std']
347 lutStd['LAMBDAB'][:] = lutCat[0]['lambdaB']
348 lutStd['ATMLAMBDA'][:] = lutCat[0]['atmLambda'][:]
349 lutStd['ATMSTDTRANS'][:] = lutCat[0]['atmStdTrans'][:]
350
351 lutTypes = [row['luttype'] for row in lutCat]
352
353 # And the flattened look-up-table
354 lutFlat = np.zeros(lutCat[0]['lut'].size, dtype=[('I0', 'f4'),
355 ('I1', 'f4')])
356
357 lutFlat['I0'][:] = lutCat[lutTypes.index('I0')]['lut'][:]
358 lutFlat['I1'][:] = lutCat[lutTypes.index('I1')]['lut'][:]
359
360 lutDerivFlat = np.zeros(lutCat[0]['lut'].size, dtype=[('D_LNPWV', 'f4'),
361 ('D_O3', 'f4'),
362 ('D_LNTAU', 'f4'),
363 ('D_ALPHA', 'f4'),
364 ('D_SECZENITH', 'f4'),
365 ('D_LNPWV_I1', 'f4'),
366 ('D_O3_I1', 'f4'),
367 ('D_LNTAU_I1', 'f4'),
368 ('D_ALPHA_I1', 'f4'),
369 ('D_SECZENITH_I1', 'f4')])
370
371 for name in lutDerivFlat.dtype.names:
372 lutDerivFlat[name][:] = lutCat[lutTypes.index(name)]['lut'][:]
373
374 # The fgcm.FgcmLUT() class copies all the LUT information into special
375 # shared memory objects that will not blow up the memory usage when used
376 # with python multiprocessing. Once all the numbers are copied, the
377 # references to the temporary objects (lutCat, lutFlat, lutDerivFlat)
378 # will fall out of scope and can be cleaned up by the garbage collector.
379 fgcmLut = fgcm.FgcmLUT(lutIndexVals, lutFlat, lutDerivFlat, lutStd,
380 filterToBand=physicalFilterMap)
381
382 return fgcmLut, lutIndexVals, lutStd
383
384
386 """
387 Translate the FGCM visit catalog to an fgcm-compatible object
388
389 Parameters
390 ----------
391 visitCat: `lsst.afw.table.BaseCatalog`
392 FGCM visitCat from `lsst.fgcmcal.FgcmBuildStarsTask`
393
394 Returns
395 -------
396 fgcmExpInfo: `numpy.ndarray`
397 Numpy array for visit information for FGCM
398
399 Notes
400 -----
401 After running this code, it is wise to `del visitCat` to clear the memory.
402 """
403
404 fgcmExpInfo = np.zeros(len(visitCat), dtype=[('VISIT', 'i8'),
405 ('MJD', 'f8'),
406 ('EXPTIME', 'f8'),
407 ('PSFSIGMA', 'f8'),
408 ('DELTA_APER', 'f8'),
409 ('SKYBACKGROUND', 'f8'),
410 ('DEEPFLAG', 'i2'),
411 ('TELHA', 'f8'),
412 ('TELRA', 'f8'),
413 ('TELDEC', 'f8'),
414 ('TELROT', 'f8'),
415 ('PMB', 'f8'),
416 ('FILTERNAME', 'a50')])
417 fgcmExpInfo['VISIT'][:] = visitCat['visit']
418 fgcmExpInfo['MJD'][:] = visitCat['mjd']
419 fgcmExpInfo['EXPTIME'][:] = visitCat['exptime']
420 fgcmExpInfo['DEEPFLAG'][:] = visitCat['deepFlag']
421 fgcmExpInfo['TELHA'][:] = visitCat['telha']
422 fgcmExpInfo['TELRA'][:] = visitCat['telra']
423 fgcmExpInfo['TELDEC'][:] = visitCat['teldec']
424 fgcmExpInfo['TELROT'][:] = visitCat['telrot']
425 fgcmExpInfo['PMB'][:] = visitCat['pmb']
426 fgcmExpInfo['PSFSIGMA'][:] = visitCat['psfSigma']
427 fgcmExpInfo['DELTA_APER'][:] = visitCat['deltaAper']
428 fgcmExpInfo['SKYBACKGROUND'][:] = visitCat['skyBackground']
429 # Note that we have to go through asAstropy() to get a string
430 # array out of an afwTable. This is faster than a row-by-row loop.
431 fgcmExpInfo['FILTERNAME'][:] = visitCat.asAstropy()['physicalFilter']
432
433 return fgcmExpInfo
434
435
437 """
438 Compute the median pixel scale in the camera
439
440 Returns
441 -------
442 pixelScale: `float`
443 Average pixel scale (arcsecond) over the camera
444 """
445
446 boresight = geom.SpherePoint(180.0*geom.degrees, 0.0*geom.degrees)
447 orientation = 0.0*geom.degrees
448 flipX = False
449
450 # Create a temporary visitInfo for input to createInitialSkyWcs
451 visitInfo = afwImage.VisitInfo(boresightRaDec=boresight,
452 boresightRotAngle=orientation,
453 rotType=afwImage.RotType.SKY)
454
455 pixelScales = np.zeros(len(camera))
456 for i, detector in enumerate(camera):
457 wcs = createInitialSkyWcs(visitInfo, detector, flipX)
458 pixelScales[i] = wcs.getPixelScale(detector.getBBox().getCenter()).asArcseconds()
459
460 ok, = np.where(pixelScales > 0.0)
461 return np.median(pixelScales[ok])
462
463
465 """
466 Compute the approximate pixel area bounded fields from the camera
467 geometry.
468
469 Parameters
470 ----------
471 camera: `lsst.afw.cameraGeom.Camera`
472
473 Returns
474 -------
475 approxPixelAreaFields: `dict`
476 Dictionary of approximate area fields, keyed with detector ID
477 """
478
479 areaScaling = 1. / computeReferencePixelScale(camera)**2.
480
481 # Generate fake WCSs centered at 180/0 to avoid the RA=0/360 problem,
482 # since we are looking for relative scales
483 boresight = geom.SpherePoint(180.0*geom.degrees, 0.0*geom.degrees)
484
485 flipX = False
486 # Create a temporary visitInfo for input to createInitialSkyWcs
487 # The orientation does not matter for the area computation
488 visitInfo = afwImage.VisitInfo(boresightRaDec=boresight,
489 boresightRotAngle=0.0*geom.degrees,
490 rotType=afwImage.RotType.SKY)
491
492 approxPixelAreaFields = {}
493
494 for i, detector in enumerate(camera):
495 key = detector.getId()
496
497 wcs = createInitialSkyWcs(visitInfo, detector, flipX)
498 bbox = detector.getBBox()
499
500 areaField = afwMath.PixelAreaBoundedField(bbox, wcs,
501 unit=geom.arcseconds, scaling=areaScaling)
502 approxAreaField = afwMath.ChebyshevBoundedField.approximate(areaField)
503
504 approxPixelAreaFields[key] = approxAreaField
505
506 return approxPixelAreaFields
507
508
509def makeZptSchema(superStarChebyshevSize, zptChebyshevSize):
510 """
511 Make the zeropoint schema
512
513 Parameters
514 ----------
515 superStarChebyshevSize: `int`
516 Length of the superstar chebyshev array
517 zptChebyshevSize: `int`
518 Length of the zeropoint chebyshev array
519
520 Returns
521 -------
522 zptSchema: `lsst.afw.table.schema`
523 """
524
525 zptSchema = afwTable.Schema()
526
527 zptSchema.addField('visit', type=np.int64, doc='Visit number')
528 zptSchema.addField('detector', type=np.int32, doc='Detector ID number')
529 zptSchema.addField('fgcmFlag', type=np.int32, doc=('FGCM flag value: '
530 '1: Photometric, used in fit; '
531 '2: Photometric, not used in fit; '
532 '4: Non-photometric, on partly photometric night; '
533 '8: Non-photometric, on non-photometric night; '
534 '16: No zeropoint could be determined; '
535 '32: Too few stars for reliable gray computation'))
536 zptSchema.addField('fgcmZpt', type=np.float64, doc='FGCM zeropoint (center of CCD)')
537 zptSchema.addField('fgcmZptErr', type=np.float64,
538 doc='Error on zeropoint, estimated from repeatability + number of obs')
539 zptSchema.addField('fgcmfZptChebXyMax', type='ArrayD', size=2,
540 doc='maximum x/maximum y to scale to apply chebyshev parameters')
541 zptSchema.addField('fgcmfZptCheb', type='ArrayD',
542 size=zptChebyshevSize,
543 doc='Chebyshev parameters (flattened) for zeropoint')
544 zptSchema.addField('fgcmfZptSstarCheb', type='ArrayD',
545 size=superStarChebyshevSize,
546 doc='Chebyshev parameters (flattened) for superStarFlat')
547 zptSchema.addField('fgcmI0', type=np.float64, doc='Integral of the passband')
548 zptSchema.addField('fgcmI10', type=np.float64, doc='Normalized chromatic integral')
549 zptSchema.addField('fgcmR0', type=np.float64,
550 doc='Retrieved i0 integral, estimated from stars (only for flag 1)')
551 zptSchema.addField('fgcmR10', type=np.float64,
552 doc='Retrieved i10 integral, estimated from stars (only for flag 1)')
553 zptSchema.addField('fgcmGry', type=np.float64,
554 doc='Estimated gray extinction relative to atmospheric solution; '
555 'only for fgcmFlag <= 4 (see fgcmFlag) ')
556 zptSchema.addField('fgcmDeltaChrom', type=np.float64,
557 doc='Mean chromatic correction for stars in this ccd; '
558 'only for fgcmFlag <= 4 (see fgcmFlag)')
559 zptSchema.addField('fgcmZptVar', type=np.float64, doc='Variance of zeropoint over ccd')
560 zptSchema.addField('fgcmTilings', type=np.float64,
561 doc='Number of photometric tilings used for solution for ccd')
562 zptSchema.addField('fgcmFpGry', type=np.float64,
563 doc='Average gray extinction over the full focal plane '
564 '(same for all ccds in a visit)')
565 zptSchema.addField('fgcmFpGryBlue', type=np.float64,
566 doc='Average gray extinction over the full focal plane '
567 'for 25% bluest stars')
568 zptSchema.addField('fgcmFpGryBlueErr', type=np.float64,
569 doc='Error on Average gray extinction over the full focal plane '
570 'for 25% bluest stars')
571 zptSchema.addField('fgcmFpGryRed', type=np.float64,
572 doc='Average gray extinction over the full focal plane '
573 'for 25% reddest stars')
574 zptSchema.addField('fgcmFpGryRedErr', type=np.float64,
575 doc='Error on Average gray extinction over the full focal plane '
576 'for 25% reddest stars')
577 zptSchema.addField('fgcmFpVar', type=np.float64,
578 doc='Variance of gray extinction over the full focal plane '
579 '(same for all ccds in a visit)')
580 zptSchema.addField('fgcmDust', type=np.float64,
581 doc='Gray dust extinction from the primary/corrector'
582 'at the time of the exposure')
583 zptSchema.addField('fgcmFlat', type=np.float64, doc='Superstarflat illumination correction')
584 zptSchema.addField('fgcmAperCorr', type=np.float64, doc='Aperture correction estimated by fgcm')
585 zptSchema.addField('fgcmDeltaMagBkg', type=np.float64,
586 doc=('Local background correction from brightest percentile '
587 '(value set by deltaMagBkgOffsetPercentile) calibration '
588 'stars.'))
589 zptSchema.addField('exptime', type=np.float32, doc='Exposure time')
590 zptSchema.addField('filtername', type=str, size=30, doc='Filter name')
591
592 return zptSchema
593
594
595def makeZptCat(zptSchema, zpStruct):
596 """
597 Make the zeropoint catalog for persistence
598
599 Parameters
600 ----------
601 zptSchema: `lsst.afw.table.Schema`
602 Zeropoint catalog schema
603 zpStruct: `numpy.ndarray`
604 Zeropoint structure from fgcm
605
606 Returns
607 -------
608 zptCat: `afwTable.BaseCatalog`
609 Zeropoint catalog for persistence
610 """
611
612 zptCat = afwTable.BaseCatalog(zptSchema)
613 zptCat.reserve(zpStruct.size)
614
615 for filterName in zpStruct['FILTERNAME']:
616 rec = zptCat.addNew()
617 rec['filtername'] = filterName.decode('utf-8')
618
619 zptCat['visit'][:] = zpStruct[FGCM_EXP_FIELD]
620 zptCat['detector'][:] = zpStruct[FGCM_CCD_FIELD]
621 zptCat['fgcmFlag'][:] = zpStruct['FGCM_FLAG']
622 zptCat['fgcmZpt'][:] = zpStruct['FGCM_ZPT']
623 zptCat['fgcmZptErr'][:] = zpStruct['FGCM_ZPTERR']
624 zptCat['fgcmfZptChebXyMax'][:, :] = zpStruct['FGCM_FZPT_XYMAX']
625 zptCat['fgcmfZptCheb'][:, :] = zpStruct['FGCM_FZPT_CHEB']
626 zptCat['fgcmfZptSstarCheb'][:, :] = zpStruct['FGCM_FZPT_SSTAR_CHEB']
627 zptCat['fgcmI0'][:] = zpStruct['FGCM_I0']
628 zptCat['fgcmI10'][:] = zpStruct['FGCM_I10']
629 zptCat['fgcmR0'][:] = zpStruct['FGCM_R0']
630 zptCat['fgcmR10'][:] = zpStruct['FGCM_R10']
631 zptCat['fgcmGry'][:] = zpStruct['FGCM_GRY']
632 zptCat['fgcmDeltaChrom'][:] = zpStruct['FGCM_DELTACHROM']
633 zptCat['fgcmZptVar'][:] = zpStruct['FGCM_ZPTVAR']
634 zptCat['fgcmTilings'][:] = zpStruct['FGCM_TILINGS']
635 zptCat['fgcmFpGry'][:] = zpStruct['FGCM_FPGRY']
636 zptCat['fgcmFpGryBlue'][:] = zpStruct['FGCM_FPGRY_CSPLIT'][:, 0]
637 zptCat['fgcmFpGryBlueErr'][:] = zpStruct['FGCM_FPGRY_CSPLITERR'][:, 0]
638 zptCat['fgcmFpGryRed'][:] = zpStruct['FGCM_FPGRY_CSPLIT'][:, 2]
639 zptCat['fgcmFpGryRedErr'][:] = zpStruct['FGCM_FPGRY_CSPLITERR'][:, 2]
640 zptCat['fgcmFpVar'][:] = zpStruct['FGCM_FPVAR']
641 zptCat['fgcmDust'][:] = zpStruct['FGCM_DUST']
642 zptCat['fgcmFlat'][:] = zpStruct['FGCM_FLAT']
643 zptCat['fgcmAperCorr'][:] = zpStruct['FGCM_APERCORR']
644 zptCat['fgcmDeltaMagBkg'][:] = zpStruct['FGCM_DELTAMAGBKG']
645 zptCat['exptime'][:] = zpStruct['EXPTIME']
646
647 return zptCat
648
649
651 """
652 Make the atmosphere schema
653
654 Returns
655 -------
656 atmSchema: `lsst.afw.table.Schema`
657 """
658
659 atmSchema = afwTable.Schema()
660
661 atmSchema.addField('visit', type=np.int64, doc='Visit number')
662 atmSchema.addField('pmb', type=np.float64, doc='Barometric pressure (mb)')
663 atmSchema.addField('pwv', type=np.float64, doc='Water vapor (mm)')
664 atmSchema.addField('tau', type=np.float64, doc='Aerosol optical depth')
665 atmSchema.addField('alpha', type=np.float64, doc='Aerosol slope')
666 atmSchema.addField('o3', type=np.float64, doc='Ozone (dobson)')
667 atmSchema.addField('secZenith', type=np.float64, doc='Secant(zenith) (~ airmass)')
668 atmSchema.addField('cTrans', type=np.float64, doc='Transmission correction factor')
669 atmSchema.addField('lamStd', type=np.float64, doc='Wavelength for transmission correction')
670
671 return atmSchema
672
673
674def makeAtmCat(atmSchema, atmStruct):
675 """
676 Make the atmosphere catalog for persistence
677
678 Parameters
679 ----------
680 atmSchema: `lsst.afw.table.Schema`
681 Atmosphere catalog schema
682 atmStruct: `numpy.ndarray`
683 Atmosphere structure from fgcm
684
685 Returns
686 -------
687 atmCat: `lsst.afw.table.BaseCatalog`
688 Atmosphere catalog for persistence
689 """
690
691 atmCat = afwTable.BaseCatalog(atmSchema)
692 atmCat.resize(atmStruct.size)
693
694 atmCat['visit'][:] = atmStruct['VISIT']
695 atmCat['pmb'][:] = atmStruct['PMB']
696 atmCat['pwv'][:] = atmStruct['PWV']
697 atmCat['tau'][:] = atmStruct['TAU']
698 atmCat['alpha'][:] = atmStruct['ALPHA']
699 atmCat['o3'][:] = atmStruct['O3']
700 atmCat['secZenith'][:] = atmStruct['SECZENITH']
701 atmCat['cTrans'][:] = atmStruct['CTRANS']
702 atmCat['lamStd'][:] = atmStruct['LAMSTD']
703
704 return atmCat
705
706
707def makeStdSchema(nBands):
708 """
709 Make the standard star schema
710
711 Parameters
712 ----------
713 nBands: `int`
714 Number of bands in standard star catalog
715
716 Returns
717 -------
718 stdSchema: `lsst.afw.table.Schema`
719 """
720
721 stdSchema = afwTable.SimpleTable.makeMinimalSchema()
722 stdSchema.addField('ngood', type='ArrayI', doc='Number of good observations',
723 size=nBands)
724 stdSchema.addField('ntotal', type='ArrayI', doc='Number of total observations',
725 size=nBands)
726 stdSchema.addField('mag_std_noabs', type='ArrayF',
727 doc='Standard magnitude (no absolute calibration)',
728 size=nBands)
729 stdSchema.addField('magErr_std', type='ArrayF',
730 doc='Standard magnitude error',
731 size=nBands)
732 stdSchema.addField('npsfcand', type='ArrayI',
733 doc='Number of observations flagged as psf candidates',
734 size=nBands)
735 stdSchema.addField('delta_aper', type='ArrayF',
736 doc='Delta mag (small - large aperture)',
737 size=nBands)
738
739 return stdSchema
740
741
742def makeStdCat(stdSchema, stdStruct, goodBands):
743 """
744 Make the standard star catalog for persistence
745
746 Parameters
747 ----------
748 stdSchema: `lsst.afw.table.Schema`
749 Standard star catalog schema
750 stdStruct: `numpy.ndarray`
751 Standard star structure in FGCM format
752 goodBands: `list`
753 List of good band names used in stdStruct
754
755 Returns
756 -------
757 stdCat: `lsst.afw.table.BaseCatalog`
758 Standard star catalog for persistence
759 """
760
761 stdCat = afwTable.SimpleCatalog(stdSchema)
762 stdCat.resize(stdStruct.size)
763
764 stdCat['id'][:] = stdStruct['FGCM_ID']
765 stdCat['coord_ra'][:] = stdStruct['RA'] * geom.degrees
766 stdCat['coord_dec'][:] = stdStruct['DEC'] * geom.degrees
767 stdCat['ngood'][:, :] = stdStruct['NGOOD'][:, :]
768 stdCat['ntotal'][:, :] = stdStruct['NTOTAL'][:, :]
769 stdCat['mag_std_noabs'][:, :] = stdStruct['MAG_STD'][:, :]
770 stdCat['magErr_std'][:, :] = stdStruct['MAGERR_STD'][:, :]
771 if 'NPSFCAND' in stdStruct.dtype.names:
772 stdCat['npsfcand'][:, :] = stdStruct['NPSFCAND'][:, :]
773 stdCat['delta_aper'][:, :] = stdStruct['DELTA_APER'][:, :]
774
775 md = PropertyList()
776 md.set("BANDS", list(goodBands))
777 stdCat.setMetadata(md)
778
779 return stdCat
780
781
783 """
784 Compute the radius associated with a CircularApertureFlux or ApFlux field.
785
786 Parameters
787 ----------
788 fluxField : `str`
789 CircularApertureFlux or ApFlux
790
791 Returns
792 -------
793 apertureRadius : `float`
794 Radius of the aperture field, in pixels.
795
796 Raises
797 ------
798 RuntimeError: Raised if flux field is not a CircularApertureFlux,
799 ApFlux, or apFlux.
800 """
801 # TODO: Move this method to more general stack method in DM-25775
802 m = re.search(r'(CircularApertureFlux|ApFlux|apFlux)_(\d+)_(\d+)_', fluxField)
803
804 if m is None:
805 raise RuntimeError(f"Flux field {fluxField} does not correspond to a CircularApertureFlux or ApFlux")
806
807 apertureRadius = float(m.groups()[1]) + float(m.groups()[2])/10.
808
809 return apertureRadius
810
811
812def extractReferenceMags(refStars, bands, filterMap):
813 """
814 Extract reference magnitudes from refStars for given bands and
815 associated filterMap.
816
817 Parameters
818 ----------
819 refStars : `astropy.table.Table` or `lsst.afw.table.BaseCatalog`
820 FGCM reference star catalog.
821 bands : `list`
822 List of bands for calibration.
823 filterMap: `dict`
824 FGCM mapping of filter to band.
825
826 Returns
827 -------
828 refMag : `np.ndarray`
829 nstar x nband array of reference magnitudes.
830 refMagErr : `np.ndarray`
831 nstar x nband array of reference magnitude errors.
832 """
833 hasAstropyMeta = False
834 try:
835 meta = refStars.meta
836 hasAstropyMeta = True
837 except AttributeError:
838 meta = refStars.getMetadata()
839
840 if 'FILTERNAMES' in meta:
841 if hasAstropyMeta:
842 filternames = meta['FILTERNAMES']
843 else:
844 filternames = meta.getArray('FILTERNAMES')
845
846 # The reference catalog that fgcm wants has one entry per band
847 # in the config file
848 refMag = np.zeros((len(refStars), len(bands)),
849 dtype=refStars['refMag'].dtype) + 99.0
850 refMagErr = np.zeros_like(refMag) + 99.0
851 for i, filtername in enumerate(filternames):
852 # We are allowed to run the fit configured so that we do not
853 # use every column in the reference catalog.
854 try:
855 band = filterMap[filtername]
856 except KeyError:
857 continue
858 try:
859 ind = bands.index(band)
860 except ValueError:
861 continue
862
863 refMag[:, ind] = refStars['refMag'][:, i]
864 refMagErr[:, ind] = refStars['refMagErr'][:, i]
865 else:
866 raise RuntimeError("FGCM reference stars missing FILTERNAMES metadata.")
867
868 return refMag, refMagErr
869
870
871def lookupStaticCalibrations(datasetType, registry, quantumDataId, collections):
872 # For static calibrations, we search with a timespan that has unbounded
873 # begin and end; we'll get an error if there's more than one match (because
874 # then it's not static).
875 timespan = Timespan(begin=None, end=None)
876 result = []
877 # First iterate over all of the data IDs for this dataset type that are
878 # consistent with the quantum data ID.
879 for dataId in registry.queryDataIds(datasetType.dimensions, dataId=quantumDataId):
880 # Find the dataset with this data ID using the unbounded timespan.
881 if ref := registry.findDataset(datasetType, dataId, collections=collections, timespan=timespan):
882 result.append(ref)
883 return result
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
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
Class for storing ordered metadata with comments.
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
translateVisitCatalog(visitCat)
Definition utilities.py:385
computeApertureRadiusFromName(fluxField)
Definition utilities.py:782
makeStdCat(stdSchema, stdStruct, goodBands)
Definition utilities.py:742
makeZptCat(zptSchema, zpStruct)
Definition utilities.py:595
lookupStaticCalibrations(datasetType, registry, quantumDataId, collections)
Definition utilities.py:871
makeConfigDict(config, log, camera, maxIter, resetFitParameters, outputZeropoints, lutFilterNames, tract=None, nCore=1, doPlots=False)
Definition utilities.py:49
computeReferencePixelScale(camera)
Definition utilities.py:436
makeAtmCat(atmSchema, atmStruct)
Definition utilities.py:674
translateFgcmLut(lutCat, physicalFilterMap)
Definition utilities.py:255
computeApproxPixelAreaFields(camera)
Definition utilities.py:464
makeZptSchema(superStarChebyshevSize, zptChebyshevSize)
Definition utilities.py:509
extractReferenceMags(refStars, bands, filterMap)
Definition utilities.py:812