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