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