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