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
fgcmOutputProducts.py
Go to the documentation of this file.
1 # See COPYRIGHT file at the top of the source tree.
2 #
3 # This file is part of fgcmcal.
4 #
5 # Developed for the LSST Data Management System.
6 # This product includes software developed by the LSST Project
7 # (https://www.lsst.org).
8 # See the COPYRIGHT file at the top-level directory of this distribution
9 # for details of code ownership.
10 #
11 # This program is free software: you can redistribute it and/or modify
12 # it under the terms of the GNU General Public License as published by
13 # the Free Software Foundation, either version 3 of the License, or
14 # (at your option) any later version.
15 #
16 # This program is distributed in the hope that it will be useful,
17 # but WITHOUT ANY WARRANTY; without even the implied warranty of
18 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 # GNU General Public License for more details.
20 #
21 # You should have received a copy of the GNU General Public License
22 # along with this program. If not, see <https://www.gnu.org/licenses/>.
23 """Make the final fgcmcal output products.
24 
25 This task takes the final output from fgcmFitCycle and produces the following
26 outputs for use in the DM stack: the FGCM standard stars in a reference
27 catalog format; the model atmospheres in "transmission_atmosphere_fgcm"
28 format; and the zeropoints in "fgcm_photoCalib" format. Optionally, the
29 task can transfer the 'absolute' calibration from a reference catalog
30 to put the fgcm standard stars in units of Jansky. This is accomplished
31 by matching stars in a sample of healpix pixels, and applying the median
32 offset per band.
33 """
34 import sys
35 import traceback
36 import copy
37 
38 import numpy as np
39 import healpy as hp
40 import esutil
41 from astropy import units
42 
43 import lsst.daf.base as dafBase
44 import lsst.pex.config as pexConfig
45 import lsst.pipe.base as pipeBase
46 from lsst.pipe.base import connectionTypes
47 from lsst.afw.image import TransmissionCurve
48 from lsst.meas.algorithms import LoadIndexedReferenceObjectsTask
49 from lsst.meas.algorithms import ReferenceObjectLoader
50 from lsst.pipe.tasks.photoCal import PhotoCalTask
51 import lsst.geom
52 import lsst.afw.image as afwImage
53 import lsst.afw.math as afwMath
54 import lsst.afw.table as afwTable
55 from lsst.meas.algorithms import IndexerRegistry
56 from lsst.meas.algorithms import DatasetConfig
57 from lsst.meas.algorithms.ingestIndexReferenceTask import addRefCatMetadata
58 from lsst.utils.timer import timeMethod
59 
60 from .utilities import computeApproxPixelAreaFields
61 from .utilities import lookupStaticCalibrations
62 from .utilities import FGCM_ILLEGAL_VALUE
63 
64 import fgcm
65 
66 __all__ = ['FgcmOutputProductsConfig', 'FgcmOutputProductsTask', 'FgcmOutputProductsRunner']
67 
68 
69 class FgcmOutputProductsConnections(pipeBase.PipelineTaskConnections,
70  dimensions=("instrument",),
71  defaultTemplates={"cycleNumber": "0"}):
72  camera = connectionTypes.PrerequisiteInput(
73  doc="Camera instrument",
74  name="camera",
75  storageClass="Camera",
76  dimensions=("instrument",),
77  lookupFunction=lookupStaticCalibrations,
78  isCalibration=True,
79  )
80 
81  fgcmLookUpTable = connectionTypes.PrerequisiteInput(
82  doc=("Atmosphere + instrument look-up-table for FGCM throughput and "
83  "chromatic corrections."),
84  name="fgcmLookUpTable",
85  storageClass="Catalog",
86  dimensions=("instrument",),
87  deferLoad=True,
88  )
89 
90  fgcmVisitCatalog = connectionTypes.Input(
91  doc="Catalog of visit information for fgcm",
92  name="fgcmVisitCatalog",
93  storageClass="Catalog",
94  dimensions=("instrument",),
95  deferLoad=True,
96  )
97 
98  fgcmStandardStars = connectionTypes.Input(
99  doc="Catalog of standard star data from fgcm fit",
100  name="fgcmStandardStars{cycleNumber}",
101  storageClass="SimpleCatalog",
102  dimensions=("instrument",),
103  deferLoad=True,
104  )
105 
106  fgcmZeropoints = connectionTypes.Input(
107  doc="Catalog of zeropoints from fgcm fit",
108  name="fgcmZeropoints{cycleNumber}",
109  storageClass="Catalog",
110  dimensions=("instrument",),
111  deferLoad=True,
112  )
113 
114  fgcmAtmosphereParameters = connectionTypes.Input(
115  doc="Catalog of atmosphere parameters from fgcm fit",
116  name="fgcmAtmosphereParameters{cycleNumber}",
117  storageClass="Catalog",
118  dimensions=("instrument",),
119  deferLoad=True,
120  )
121 
122  refCat = connectionTypes.PrerequisiteInput(
123  doc="Reference catalog to use for photometric calibration",
124  name="cal_ref_cat",
125  storageClass="SimpleCatalog",
126  dimensions=("skypix",),
127  deferLoad=True,
128  multiple=True,
129  )
130 
131  fgcmPhotoCalib = connectionTypes.Output(
132  doc=("Per-visit photometric calibrations derived from fgcm calibration. "
133  "These catalogs use detector id for the id and are sorted for "
134  "fast lookups of a detector."),
135  name="fgcmPhotoCalibCatalog",
136  storageClass="ExposureCatalog",
137  dimensions=("instrument", "visit",),
138  multiple=True,
139  )
140 
141  fgcmTransmissionAtmosphere = connectionTypes.Output(
142  doc="Per-visit atmosphere transmission files produced from fgcm calibration",
143  name="transmission_atmosphere_fgcm",
144  storageClass="TransmissionCurve",
145  dimensions=("instrument",
146  "visit",),
147  multiple=True,
148  )
149 
150  fgcmOffsets = connectionTypes.Output(
151  doc="Per-band offsets computed from doReferenceCalibration",
152  name="fgcmReferenceCalibrationOffsets",
153  storageClass="Catalog",
154  dimensions=("instrument",),
155  multiple=False,
156  )
157 
158  def __init__(self, *, config=None):
159  super().__init__(config=config)
160 
161  if str(int(config.connections.cycleNumber)) != config.connections.cycleNumber:
162  raise ValueError("cycleNumber must be of integer format")
163  if config.connections.refCat != config.refObjLoader.ref_dataset_name:
164  raise ValueError("connections.refCat must be the same as refObjLoader.ref_dataset_name")
165 
166  if config.doRefcatOutput:
167  raise ValueError("FgcmOutputProductsTask (Gen3) does not support doRefcatOutput")
168 
169  if not config.doReferenceCalibration:
170  self.prerequisiteInputs.remove("refCat")
171  if not config.doAtmosphereOutput:
172  self.inputs.remove("fgcmAtmosphereParameters")
173  if not config.doZeropointOutput:
174  self.inputs.remove("fgcmZeropoints")
175  if not config.doReferenceCalibration:
176  self.outputs.remove("fgcmOffsets")
177 
178 
179 class FgcmOutputProductsConfig(pipeBase.PipelineTaskConfig,
180  pipelineConnections=FgcmOutputProductsConnections):
181  """Config for FgcmOutputProductsTask"""
182 
183  cycleNumber = pexConfig.Field(
184  doc="Final fit cycle from FGCM fit",
185  dtype=int,
186  default=None,
187  )
188  physicalFilterMap = pexConfig.DictField(
189  doc="Mapping from 'physicalFilter' to band.",
190  keytype=str,
191  itemtype=str,
192  default={},
193  )
194  # The following fields refer to calibrating from a reference
195  # catalog, but in the future this might need to be expanded
196  doReferenceCalibration = pexConfig.Field(
197  doc=("Transfer 'absolute' calibration from reference catalog? "
198  "This afterburner step is unnecessary if reference stars "
199  "were used in the full fit in FgcmFitCycleTask."),
200  dtype=bool,
201  default=False,
202  )
203  doRefcatOutput = pexConfig.Field(
204  doc="Output standard stars in reference catalog format",
205  dtype=bool,
206  default=True,
207  )
208  doAtmosphereOutput = pexConfig.Field(
209  doc="Output atmospheres in transmission_atmosphere_fgcm format",
210  dtype=bool,
211  default=True,
212  )
213  doZeropointOutput = pexConfig.Field(
214  doc="Output zeropoints in fgcm_photoCalib format",
215  dtype=bool,
216  default=True,
217  )
218  doComposeWcsJacobian = pexConfig.Field(
219  doc="Compose Jacobian of WCS with fgcm calibration for output photoCalib?",
220  dtype=bool,
221  default=True,
222  )
223  doApplyMeanChromaticCorrection = pexConfig.Field(
224  doc="Apply the mean chromatic correction to the zeropoints?",
225  dtype=bool,
226  default=True,
227  )
228  refObjLoader = pexConfig.ConfigurableField(
229  target=LoadIndexedReferenceObjectsTask,
230  doc="reference object loader for 'absolute' photometric calibration",
231  )
232  photoCal = pexConfig.ConfigurableField(
233  target=PhotoCalTask,
234  doc="task to perform 'absolute' calibration",
235  )
236  referencePixelizationNside = pexConfig.Field(
237  doc="Healpix nside to pixelize catalog to compare to reference catalog",
238  dtype=int,
239  default=64,
240  )
241  referencePixelizationMinStars = pexConfig.Field(
242  doc=("Minimum number of stars per healpix pixel to select for comparison"
243  "to the specified reference catalog"),
244  dtype=int,
245  default=200,
246  )
247  referenceMinMatch = pexConfig.Field(
248  doc="Minimum number of stars matched to reference catalog to be used in statistics",
249  dtype=int,
250  default=50,
251  )
252  referencePixelizationNPixels = pexConfig.Field(
253  doc=("Number of healpix pixels to sample to do comparison. "
254  "Doing too many will take a long time and not yield any more "
255  "precise results because the final number is the median offset "
256  "(per band) from the set of pixels."),
257  dtype=int,
258  default=100,
259  )
260  datasetConfig = pexConfig.ConfigField(
261  dtype=DatasetConfig,
262  doc="Configuration for writing/reading ingested catalog",
263  )
264 
265  def setDefaults(self):
266  pexConfig.Config.setDefaults(self)
267 
268  # In order to transfer the "absolute" calibration from a reference
269  # catalog to the relatively calibrated FGCM standard stars (one number
270  # per band), we use the PhotoCalTask to match stars in a sample of healpix
271  # pixels. These basic settings ensure that only well-measured, good stars
272  # from the source and reference catalogs are used for the matching.
273 
274  # applyColorTerms needs to be False if doReferenceCalibration is False,
275  # as is the new default after DM-16702
276  self.photoCal.applyColorTerms = False
277  self.photoCal.fluxField = 'instFlux'
278  self.photoCal.magErrFloor = 0.003
279  self.photoCal.match.referenceSelection.doSignalToNoise = True
280  self.photoCal.match.referenceSelection.signalToNoise.minimum = 10.0
281  self.photoCal.match.sourceSelection.doSignalToNoise = True
282  self.photoCal.match.sourceSelection.signalToNoise.minimum = 10.0
283  self.photoCal.match.sourceSelection.signalToNoise.fluxField = 'instFlux'
284  self.photoCal.match.sourceSelection.signalToNoise.errField = 'instFluxErr'
285  self.photoCal.match.sourceSelection.doFlags = True
286  self.photoCal.match.sourceSelection.flags.good = []
287  self.photoCal.match.sourceSelection.flags.bad = ['flag_badStar']
288  self.photoCal.match.sourceSelection.doUnresolved = False
289  self.datasetConfig.ref_dataset_name = 'fgcm_stars'
290  self.datasetConfig.format_version = 1
291 
292  def validate(self):
293  super().validate()
294 
295  # Force the connections to conform with cycleNumber
296  self.connections.cycleNumber = str(self.cycleNumber)
297 
298 
299 class FgcmOutputProductsRunner(pipeBase.ButlerInitializedTaskRunner):
300  """Subclass of TaskRunner for fgcmOutputProductsTask
301 
302  fgcmOutputProductsTask.run() takes one argument, the butler, and
303  does not run on any data in the repository.
304  This runner does not use any parallelization.
305  """
306 
307  @staticmethod
308  def getTargetList(parsedCmd):
309  """
310  Return a list with one element, the butler.
311  """
312  return [parsedCmd.butler]
313 
314  def __call__(self, butler):
315  """
316  Parameters
317  ----------
318  butler: `lsst.daf.persistence.Butler`
319 
320  Returns
321  -------
322  exitStatus: `list` with `pipeBase.Struct`
323  exitStatus (0: success; 1: failure)
324  if self.doReturnResults also
325  results (`np.array` with absolute zeropoint offsets)
326  """
327  task = self.TaskClass(butler=butler, config=self.config, log=self.log)
328 
329  exitStatus = 0
330  if self.doRaise:
331  results = task.runDataRef(butler)
332  else:
333  try:
334  results = task.runDataRef(butler)
335  except Exception as e:
336  exitStatus = 1
337  task.log.fatal("Failed: %s" % e)
338  if not isinstance(e, pipeBase.TaskError):
339  traceback.print_exc(file=sys.stderr)
340 
341  task.writeMetadata(butler)
342 
343  if self.doReturnResults:
344  # The results here are the zeropoint offsets for each band
345  return [pipeBase.Struct(exitStatus=exitStatus,
346  results=results)]
347  else:
348  return [pipeBase.Struct(exitStatus=exitStatus)]
349 
350  def run(self, parsedCmd):
351  """
352  Run the task, with no multiprocessing
353 
354  Parameters
355  ----------
356  parsedCmd: `lsst.pipe.base.ArgumentParser` parsed command line
357  """
358 
359  resultList = []
360 
361  if self.precall(parsedCmd):
362  targetList = self.getTargetList(parsedCmd)
363  # make sure that we only get 1
364  resultList = self(targetList[0])
365 
366  return resultList
367 
368 
369 class FgcmOutputProductsTask(pipeBase.PipelineTask, pipeBase.CmdLineTask):
370  """
371  Output products from FGCM global calibration.
372  """
373 
374  ConfigClass = FgcmOutputProductsConfig
375  RunnerClass = FgcmOutputProductsRunner
376  _DefaultName = "fgcmOutputProducts"
377 
378  def __init__(self, butler=None, **kwargs):
379  super().__init__(**kwargs)
380 
381  # no saving of metadata for now
382  def _getMetadataName(self):
383  return None
384 
385  def runQuantum(self, butlerQC, inputRefs, outputRefs):
386  dataRefDict = {}
387  dataRefDict['camera'] = butlerQC.get(inputRefs.camera)
388  dataRefDict['fgcmLookUpTable'] = butlerQC.get(inputRefs.fgcmLookUpTable)
389  dataRefDict['fgcmVisitCatalog'] = butlerQC.get(inputRefs.fgcmVisitCatalog)
390  dataRefDict['fgcmStandardStars'] = butlerQC.get(inputRefs.fgcmStandardStars)
391 
392  if self.config.doZeropointOutput:
393  dataRefDict['fgcmZeropoints'] = butlerQC.get(inputRefs.fgcmZeropoints)
394  photoCalibRefDict = {photoCalibRef.dataId.byName()['visit']:
395  photoCalibRef for photoCalibRef in outputRefs.fgcmPhotoCalib}
396 
397  if self.config.doAtmosphereOutput:
398  dataRefDict['fgcmAtmosphereParameters'] = butlerQC.get(inputRefs.fgcmAtmosphereParameters)
399  atmRefDict = {atmRef.dataId.byName()['visit']: atmRef for
400  atmRef in outputRefs.fgcmTransmissionAtmosphere}
401 
402  if self.config.doReferenceCalibration:
403  refConfig = self.config.refObjLoader
404  self.refObjLoader = ReferenceObjectLoader(dataIds=[ref.datasetRef.dataId
405  for ref in inputRefs.refCat],
406  refCats=butlerQC.get(inputRefs.refCat),
407  config=refConfig,
408  log=self.log)
409  else:
410  self.refObjLoader = None
411 
412  struct = self.run(dataRefDict, self.config.physicalFilterMap, returnCatalogs=True)
413 
414  # Output the photoCalib exposure catalogs
415  if struct.photoCalibCatalogs is not None:
416  self.log.info("Outputting photoCalib catalogs.")
417  for visit, expCatalog in struct.photoCalibCatalogs:
418  butlerQC.put(expCatalog, photoCalibRefDict[visit])
419  self.log.info("Done outputting photoCalib catalogs.")
420 
421  # Output the atmospheres
422  if struct.atmospheres is not None:
423  self.log.info("Outputting atmosphere transmission files.")
424  for visit, atm in struct.atmospheres:
425  butlerQC.put(atm, atmRefDict[visit])
426  self.log.info("Done outputting atmosphere files.")
427 
428  if self.config.doReferenceCalibration:
429  # Turn offset into simple catalog for persistence if necessary
430  schema = afwTable.Schema()
431  schema.addField('offset', type=np.float64,
432  doc="Post-process calibration offset (mag)")
433  offsetCat = afwTable.BaseCatalog(schema)
434  offsetCat.resize(len(struct.offsets))
435  offsetCat['offset'][:] = struct.offsets
436 
437  butlerQC.put(offsetCat, outputRefs.fgcmOffsets)
438 
439  return
440 
441  @timeMethod
442  def runDataRef(self, butler):
443  """
444  Make FGCM output products for use in the stack
445 
446  Parameters
447  ----------
448  butler: `lsst.daf.persistence.Butler`
449  cycleNumber: `int`
450  Final fit cycle number, override config.
451 
452  Returns
453  -------
454  offsets: `lsst.pipe.base.Struct`
455  A structure with array of zeropoint offsets
456 
457  Raises
458  ------
459  RuntimeError:
460  Raised if any one of the following is true:
461 
462  - butler cannot find "fgcmBuildStars_config" or
463  "fgcmBuildStarsTable_config".
464  - butler cannot find "fgcmFitCycle_config".
465  - "fgcmFitCycle_config" does not refer to
466  `self.config.cycleNumber`.
467  - butler cannot find "fgcmAtmosphereParameters" and
468  `self.config.doAtmosphereOutput` is `True`.
469  - butler cannot find "fgcmStandardStars" and
470  `self.config.doReferenceCalibration` is `True` or
471  `self.config.doRefcatOutput` is `True`.
472  - butler cannot find "fgcmZeropoints" and
473  `self.config.doZeropointOutput` is `True`.
474  """
475  if self.config.doReferenceCalibration:
476  # We need the ref obj loader to get the flux field
477  self.makeSubtask("refObjLoader", butler=butler)
478 
479  # Check to make sure that the fgcmBuildStars config exists, to retrieve
480  # the visit and ccd dataset tags
481  if not butler.datasetExists('fgcmBuildStarsTable_config') and \
482  not butler.datasetExists('fgcmBuildStars_config'):
483  raise RuntimeError("Cannot find fgcmBuildStarsTable_config or fgcmBuildStars_config, "
484  "which is prereq for fgcmOutputProducts")
485 
486  if butler.datasetExists('fgcmBuildStarsTable_config'):
487  fgcmBuildStarsConfig = butler.get('fgcmBuildStarsTable_config')
488  else:
489  fgcmBuildStarsConfig = butler.get('fgcmBuildStars_config')
490  visitDataRefName = fgcmBuildStarsConfig.visitDataRefName
491  ccdDataRefName = fgcmBuildStarsConfig.ccdDataRefName
492  physicalFilterMap = fgcmBuildStarsConfig.physicalFilterMap
493 
494  if self.config.doComposeWcsJacobian and not fgcmBuildStarsConfig.doApplyWcsJacobian:
495  raise RuntimeError("Cannot compose the WCS jacobian if it hasn't been applied "
496  "in fgcmBuildStarsTask.")
497 
498  if not self.config.doComposeWcsJacobian and fgcmBuildStarsConfig.doApplyWcsJacobian:
499  self.log.warning("Jacobian was applied in build-stars but doComposeWcsJacobian is not set.")
500 
501  # And make sure that the atmosphere was output properly
502  if (self.config.doAtmosphereOutput
503  and not butler.datasetExists('fgcmAtmosphereParameters', fgcmcycle=self.config.cycleNumber)):
504  raise RuntimeError(f"Atmosphere parameters are missing for cycle {self.config.cycleNumber}.")
505 
506  if not butler.datasetExists('fgcmStandardStars',
507  fgcmcycle=self.config.cycleNumber):
508  raise RuntimeError("Standard stars are missing for cycle %d." %
509  (self.config.cycleNumber))
510 
511  if (self.config.doZeropointOutput
512  and (not butler.datasetExists('fgcmZeropoints', fgcmcycle=self.config.cycleNumber))):
513  raise RuntimeError("Zeropoints are missing for cycle %d." %
514  (self.config.cycleNumber))
515 
516  dataRefDict = {}
517  # This is the _actual_ camera
518  dataRefDict['camera'] = butler.get('camera')
519  dataRefDict['fgcmLookUpTable'] = butler.dataRef('fgcmLookUpTable')
520  dataRefDict['fgcmVisitCatalog'] = butler.dataRef('fgcmVisitCatalog')
521  dataRefDict['fgcmStandardStars'] = butler.dataRef('fgcmStandardStars',
522  fgcmcycle=self.config.cycleNumber)
523 
524  if self.config.doZeropointOutput:
525  dataRefDict['fgcmZeropoints'] = butler.dataRef('fgcmZeropoints',
526  fgcmcycle=self.config.cycleNumber)
527  if self.config.doAtmosphereOutput:
528  dataRefDict['fgcmAtmosphereParameters'] = butler.dataRef('fgcmAtmosphereParameters',
529  fgcmcycle=self.config.cycleNumber)
530 
531  struct = self.run(dataRefDict, physicalFilterMap, butler=butler, returnCatalogs=False)
532 
533  if struct.photoCalibs is not None:
534  self.log.info("Outputting photoCalib files.")
535 
536  for visit, detector, physicalFilter, photoCalib in struct.photoCalibs:
537  butler.put(photoCalib, 'fgcm_photoCalib',
538  dataId={visitDataRefName: visit,
539  ccdDataRefName: detector,
540  'filter': physicalFilter})
541 
542  self.log.info("Done outputting photoCalib files.")
543 
544  if struct.atmospheres is not None:
545  self.log.info("Outputting atmosphere transmission files.")
546  for visit, atm in struct.atmospheres:
547  butler.put(atm, "transmission_atmosphere_fgcm",
548  dataId={visitDataRefName: visit})
549  self.log.info("Done outputting atmosphere transmissions.")
550 
551  return pipeBase.Struct(offsets=struct.offsets)
552 
553  def run(self, dataRefDict, physicalFilterMap, returnCatalogs=True, butler=None):
554  """Run the output products task.
555 
556  Parameters
557  ----------
558  dataRefDict : `dict`
559  All dataRefs are `lsst.daf.persistence.ButlerDataRef` (gen2) or
560  `lsst.daf.butler.DeferredDatasetHandle` (gen3)
561  dataRef dictionary with keys:
562 
563  ``"camera"``
564  Camera object (`lsst.afw.cameraGeom.Camera`)
565  ``"fgcmLookUpTable"``
566  dataRef for the FGCM look-up table.
567  ``"fgcmVisitCatalog"``
568  dataRef for visit summary catalog.
569  ``"fgcmStandardStars"``
570  dataRef for the output standard star catalog.
571  ``"fgcmZeropoints"``
572  dataRef for the zeropoint data catalog.
573  ``"fgcmAtmosphereParameters"``
574  dataRef for the atmosphere parameter catalog.
575  ``"fgcmBuildStarsTableConfig"``
576  Config for `lsst.fgcmcal.fgcmBuildStarsTableTask`.
577  physicalFilterMap : `dict`
578  Dictionary of mappings from physical filter to FGCM band.
579  returnCatalogs : `bool`, optional
580  Return photoCalibs as per-visit exposure catalogs.
581  butler : `lsst.daf.persistence.Butler`, optional
582  Gen2 butler used for reference star outputs
583 
584  Returns
585  -------
586  retStruct : `lsst.pipe.base.Struct`
587  Output structure with keys:
588 
589  offsets : `np.ndarray`
590  Final reference offsets, per band.
591  atmospheres : `generator` [(`int`, `lsst.afw.image.TransmissionCurve`)]
592  Generator that returns (visit, transmissionCurve) tuples.
593  photoCalibs : `generator` [(`int`, `int`, `str`, `lsst.afw.image.PhotoCalib`)]
594  Generator that returns (visit, ccd, filtername, photoCalib) tuples.
595  (returned if returnCatalogs is False).
596  photoCalibCatalogs : `generator` [(`int`, `lsst.afw.table.ExposureCatalog`)]
597  Generator that returns (visit, exposureCatalog) tuples.
598  (returned if returnCatalogs is True).
599  """
600  stdCat = dataRefDict['fgcmStandardStars'].get()
601  md = stdCat.getMetadata()
602  bands = md.getArray('BANDS')
603 
604  if self.config.doReferenceCalibration:
605  lutCat = dataRefDict['fgcmLookUpTable'].get()
606  offsets = self._computeReferenceOffsets(stdCat, lutCat, physicalFilterMap, bands)
607  else:
608  offsets = np.zeros(len(bands))
609 
610  # This is Gen2 only, and requires the butler.
611  if self.config.doRefcatOutput and butler is not None:
612  self._outputStandardStars(butler, stdCat, offsets, bands, self.config.datasetConfig)
613 
614  del stdCat
615 
616  if self.config.doZeropointOutput:
617  zptCat = dataRefDict['fgcmZeropoints'].get()
618  visitCat = dataRefDict['fgcmVisitCatalog'].get()
619 
620  pcgen = self._outputZeropoints(dataRefDict['camera'], zptCat, visitCat, offsets, bands,
621  physicalFilterMap, returnCatalogs=returnCatalogs)
622  else:
623  pcgen = None
624 
625  if self.config.doAtmosphereOutput:
626  atmCat = dataRefDict['fgcmAtmosphereParameters'].get()
627  atmgen = self._outputAtmospheres(dataRefDict, atmCat)
628  else:
629  atmgen = None
630 
631  retStruct = pipeBase.Struct(offsets=offsets,
632  atmospheres=atmgen)
633  if returnCatalogs:
634  retStruct.photoCalibCatalogs = pcgen
635  else:
636  retStruct.photoCalibs = pcgen
637 
638  return retStruct
639 
640  def generateTractOutputProducts(self, dataRefDict, tract,
641  visitCat, zptCat, atmCat, stdCat,
642  fgcmBuildStarsConfig,
643  returnCatalogs=True,
644  butler=None):
645  """
646  Generate the output products for a given tract, as specified in the config.
647 
648  This method is here to have an alternate entry-point for
649  FgcmCalibrateTract.
650 
651  Parameters
652  ----------
653  dataRefDict : `dict`
654  All dataRefs are `lsst.daf.persistence.ButlerDataRef` (gen2) or
655  `lsst.daf.butler.DeferredDatasetHandle` (gen3)
656  dataRef dictionary with keys:
657 
658  ``"camera"``
659  Camera object (`lsst.afw.cameraGeom.Camera`)
660  ``"fgcmLookUpTable"``
661  dataRef for the FGCM look-up table.
662  tract : `int`
663  Tract number
664  visitCat : `lsst.afw.table.BaseCatalog`
665  FGCM visitCat from `FgcmBuildStarsTask`
666  zptCat : `lsst.afw.table.BaseCatalog`
667  FGCM zeropoint catalog from `FgcmFitCycleTask`
668  atmCat : `lsst.afw.table.BaseCatalog`
669  FGCM atmosphere parameter catalog from `FgcmFitCycleTask`
670  stdCat : `lsst.afw.table.SimpleCatalog`
671  FGCM standard star catalog from `FgcmFitCycleTask`
672  fgcmBuildStarsConfig : `lsst.fgcmcal.FgcmBuildStarsConfig`
673  Configuration object from `FgcmBuildStarsTask`
674  returnCatalogs : `bool`, optional
675  Return photoCalibs as per-visit exposure catalogs.
676  butler: `lsst.daf.persistence.Butler`, optional
677  Gen2 butler used for reference star outputs
678 
679  Returns
680  -------
681  retStruct : `lsst.pipe.base.Struct`
682  Output structure with keys:
683 
684  offsets : `np.ndarray`
685  Final reference offsets, per band.
686  atmospheres : `generator` [(`int`, `lsst.afw.image.TransmissionCurve`)]
687  Generator that returns (visit, transmissionCurve) tuples.
688  photoCalibs : `generator` [(`int`, `int`, `str`, `lsst.afw.image.PhotoCalib`)]
689  Generator that returns (visit, ccd, filtername, photoCalib) tuples.
690  (returned if returnCatalogs is False).
691  photoCalibCatalogs : `generator` [(`int`, `lsst.afw.table.ExposureCatalog`)]
692  Generator that returns (visit, exposureCatalog) tuples.
693  (returned if returnCatalogs is True).
694  """
695  physicalFilterMap = fgcmBuildStarsConfig.physicalFilterMap
696 
697  md = stdCat.getMetadata()
698  bands = md.getArray('BANDS')
699 
700  if self.config.doComposeWcsJacobian and not fgcmBuildStarsConfig.doApplyWcsJacobian:
701  raise RuntimeError("Cannot compose the WCS jacobian if it hasn't been applied "
702  "in fgcmBuildStarsTask.")
703 
704  if not self.config.doComposeWcsJacobian and fgcmBuildStarsConfig.doApplyWcsJacobian:
705  self.log.warning("Jacobian was applied in build-stars but doComposeWcsJacobian is not set.")
706 
707  if self.config.doReferenceCalibration:
708  lutCat = dataRefDict['fgcmLookUpTable'].get()
709  offsets = self._computeReferenceOffsets(stdCat, lutCat, bands, physicalFilterMap)
710  else:
711  offsets = np.zeros(len(bands))
712 
713  if self.config.doRefcatOutput and butler is not None:
714  # Create a special config that has the tract number in it
715  datasetConfig = copy.copy(self.config.datasetConfig)
716  datasetConfig.ref_dataset_name = '%s_%d' % (self.config.datasetConfig.ref_dataset_name,
717  tract)
718  self._outputStandardStars(butler, stdCat, offsets, bands, datasetConfig)
719 
720  if self.config.doZeropointOutput:
721  pcgen = self._outputZeropoints(dataRefDict['camera'], zptCat, visitCat, offsets, bands,
722  physicalFilterMap, returnCatalogs=returnCatalogs)
723  else:
724  pcgen = None
725 
726  if self.config.doAtmosphereOutput:
727  atmgen = self._outputAtmospheres(dataRefDict, atmCat)
728  else:
729  atmgen = None
730 
731  retStruct = pipeBase.Struct(offsets=offsets,
732  atmospheres=atmgen)
733  if returnCatalogs:
734  retStruct.photoCalibCatalogs = pcgen
735  else:
736  retStruct.photoCalibs = pcgen
737 
738  return retStruct
739 
740  def _computeReferenceOffsets(self, stdCat, lutCat, physicalFilterMap, bands):
741  """
742  Compute offsets relative to a reference catalog.
743 
744  This method splits the star catalog into healpix pixels
745  and computes the calibration transfer for a sample of
746  these pixels to approximate the 'absolute' calibration
747  values (on for each band) to apply to transfer the
748  absolute scale.
749 
750  Parameters
751  ----------
752  stdCat : `lsst.afw.table.SimpleCatalog`
753  FGCM standard stars
754  lutCat : `lsst.afw.table.SimpleCatalog`
755  FGCM Look-up table
756  physicalFilterMap : `dict`
757  Dictionary of mappings from physical filter to FGCM band.
758  bands : `list` [`str`]
759  List of band names from FGCM output
760  Returns
761  -------
762  offsets : `numpy.array` of floats
763  Per band zeropoint offsets
764  """
765 
766  # Only use stars that are observed in all the bands that were actually used
767  # This will ensure that we use the same healpix pixels for the absolute
768  # calibration of each band.
769  minObs = stdCat['ngood'].min(axis=1)
770 
771  goodStars = (minObs >= 1)
772  stdCat = stdCat[goodStars]
773 
774  self.log.info("Found %d stars with at least 1 good observation in each band" %
775  (len(stdCat)))
776 
777  # Associate each band with the appropriate physicalFilter and make
778  # filterLabels
779  filterLabels = []
780 
781  lutPhysicalFilters = lutCat[0]['physicalFilters'].split(',')
782  lutStdPhysicalFilters = lutCat[0]['stdPhysicalFilters'].split(',')
783  physicalFilterMapBands = list(physicalFilterMap.values())
784  physicalFilterMapFilters = list(physicalFilterMap.keys())
785  for band in bands:
786  # Find a physical filter associated from the band by doing
787  # a reverse lookup on the physicalFilterMap dict
788  physicalFilterMapIndex = physicalFilterMapBands.index(band)
789  physicalFilter = physicalFilterMapFilters[physicalFilterMapIndex]
790  # Find the appropriate fgcm standard physicalFilter
791  lutPhysicalFilterIndex = lutPhysicalFilters.index(physicalFilter)
792  stdPhysicalFilter = lutStdPhysicalFilters[lutPhysicalFilterIndex]
793  filterLabels.append(afwImage.FilterLabel(band=band,
794  physical=stdPhysicalFilter))
795 
796  # We have to make a table for each pixel with flux/fluxErr
797  # This is a temporary table generated for input to the photoCal task.
798  # These fluxes are not instFlux (they are top-of-the-atmosphere approximate and
799  # have had chromatic corrections applied to get to the standard system
800  # specified by the atmosphere/instrumental parameters), nor are they
801  # in Jansky (since they don't have a proper absolute calibration: the overall
802  # zeropoint is estimated from the telescope size, etc.)
803  sourceMapper = afwTable.SchemaMapper(stdCat.schema)
804  sourceMapper.addMinimalSchema(afwTable.SimpleTable.makeMinimalSchema())
805  sourceMapper.editOutputSchema().addField('instFlux', type=np.float64,
806  doc="instrumental flux (counts)")
807  sourceMapper.editOutputSchema().addField('instFluxErr', type=np.float64,
808  doc="instrumental flux error (counts)")
809  badStarKey = sourceMapper.editOutputSchema().addField('flag_badStar',
810  type='Flag',
811  doc="bad flag")
812 
813  # Split up the stars
814  # Note that there is an assumption here that the ra/dec coords stored
815  # on-disk are in radians, and therefore that starObs['coord_ra'] /
816  # starObs['coord_dec'] return radians when used as an array of numpy float64s.
817  theta = np.pi/2. - stdCat['coord_dec']
818  phi = stdCat['coord_ra']
819 
820  ipring = hp.ang2pix(self.config.referencePixelizationNside, theta, phi)
821  h, rev = esutil.stat.histogram(ipring, rev=True)
822 
823  gdpix, = np.where(h >= self.config.referencePixelizationMinStars)
824 
825  self.log.info("Found %d pixels (nside=%d) with at least %d good stars" %
826  (gdpix.size,
827  self.config.referencePixelizationNside,
828  self.config.referencePixelizationMinStars))
829 
830  if gdpix.size < self.config.referencePixelizationNPixels:
831  self.log.warning("Found fewer good pixels (%d) than preferred in configuration (%d)" %
832  (gdpix.size, self.config.referencePixelizationNPixels))
833  else:
834  # Sample out the pixels we want to use
835  gdpix = np.random.choice(gdpix, size=self.config.referencePixelizationNPixels, replace=False)
836 
837  results = np.zeros(gdpix.size, dtype=[('hpix', 'i4'),
838  ('nstar', 'i4', len(bands)),
839  ('nmatch', 'i4', len(bands)),
840  ('zp', 'f4', len(bands)),
841  ('zpErr', 'f4', len(bands))])
842  results['hpix'] = ipring[rev[rev[gdpix]]]
843 
844  # We need a boolean index to deal with catalogs...
845  selected = np.zeros(len(stdCat), dtype=bool)
846 
847  refFluxFields = [None]*len(bands)
848 
849  for p_index, pix in enumerate(gdpix):
850  i1a = rev[rev[pix]: rev[pix + 1]]
851 
852  # the stdCat afwTable can only be indexed with boolean arrays,
853  # and not numpy index arrays (see DM-16497). This little trick
854  # converts the index array into a boolean array
855  selected[:] = False
856  selected[i1a] = True
857 
858  for b_index, filterLabel in enumerate(filterLabels):
859  struct = self._computeOffsetOneBand(sourceMapper, badStarKey, b_index,
860  filterLabel, stdCat,
861  selected, refFluxFields)
862  results['nstar'][p_index, b_index] = len(i1a)
863  results['nmatch'][p_index, b_index] = len(struct.arrays.refMag)
864  results['zp'][p_index, b_index] = struct.zp
865  results['zpErr'][p_index, b_index] = struct.sigma
866 
867  # And compute the summary statistics
868  offsets = np.zeros(len(bands))
869 
870  for b_index, band in enumerate(bands):
871  # make configurable
872  ok, = np.where(results['nmatch'][:, b_index] >= self.config.referenceMinMatch)
873  offsets[b_index] = np.median(results['zp'][ok, b_index])
874  # use median absolute deviation to estimate Normal sigma
875  # see https://en.wikipedia.org/wiki/Median_absolute_deviation
876  madSigma = 1.4826*np.median(np.abs(results['zp'][ok, b_index] - offsets[b_index]))
877  self.log.info("Reference catalog offset for %s band: %.12f +/- %.12f",
878  band, offsets[b_index], madSigma)
879 
880  return offsets
881 
882  def _computeOffsetOneBand(self, sourceMapper, badStarKey,
883  b_index, filterLabel, stdCat, selected, refFluxFields):
884  """
885  Compute the zeropoint offset between the fgcm stdCat and the reference
886  stars for one pixel in one band
887 
888  Parameters
889  ----------
890  sourceMapper : `lsst.afw.table.SchemaMapper`
891  Mapper to go from stdCat to calibratable catalog
892  badStarKey : `lsst.afw.table.Key`
893  Key for the field with bad stars
894  b_index : `int`
895  Index of the band in the star catalog
896  filterLabel : `lsst.afw.image.FilterLabel`
897  filterLabel with band and physical filter
898  stdCat : `lsst.afw.table.SimpleCatalog`
899  FGCM standard stars
900  selected : `numpy.array(dtype=bool)`
901  Boolean array of which stars are in the pixel
902  refFluxFields : `list`
903  List of names of flux fields for reference catalog
904  """
905 
906  sourceCat = afwTable.SimpleCatalog(sourceMapper.getOutputSchema())
907  sourceCat.reserve(selected.sum())
908  sourceCat.extend(stdCat[selected], mapper=sourceMapper)
909  sourceCat['instFlux'] = 10.**(stdCat['mag_std_noabs'][selected, b_index]/(-2.5))
910  sourceCat['instFluxErr'] = (np.log(10.)/2.5)*(stdCat['magErr_std'][selected, b_index]
911  * sourceCat['instFlux'])
912  # Make sure we only use stars that have valid measurements
913  # (This is perhaps redundant with requirements above that the
914  # stars be observed in all bands, but it can't hurt)
915  badStar = (stdCat['mag_std_noabs'][selected, b_index] > 90.0)
916  for rec in sourceCat[badStar]:
917  rec.set(badStarKey, True)
918 
919  exposure = afwImage.ExposureF()
920  exposure.setFilterLabel(filterLabel)
921 
922  if refFluxFields[b_index] is None:
923  # Need to find the flux field in the reference catalog
924  # to work around limitations of DirectMatch in PhotoCal
925  ctr = stdCat[0].getCoord()
926  rad = 0.05*lsst.geom.degrees
927  refDataTest = self.refObjLoader.loadSkyCircle(ctr, rad, filterLabel.bandLabel)
928  refFluxFields[b_index] = refDataTest.fluxField
929 
930  # Make a copy of the config so that we can modify it
931  calConfig = copy.copy(self.config.photoCal.value)
932  calConfig.match.referenceSelection.signalToNoise.fluxField = refFluxFields[b_index]
933  calConfig.match.referenceSelection.signalToNoise.errField = refFluxFields[b_index] + 'Err'
934  calTask = self.config.photoCal.target(refObjLoader=self.refObjLoader,
935  config=calConfig,
936  schema=sourceCat.getSchema())
937 
938  struct = calTask.run(exposure, sourceCat)
939 
940  return struct
941 
942  def _outputStandardStars(self, butler, stdCat, offsets, bands, datasetConfig):
943  """
944  Output standard stars in indexed reference catalog format.
945  This is not currently supported in Gen3.
946 
947  Parameters
948  ----------
949  butler : `lsst.daf.persistence.Butler`
950  stdCat : `lsst.afw.table.SimpleCatalog`
951  FGCM standard star catalog from fgcmFitCycleTask
952  offsets : `numpy.array` of floats
953  Per band zeropoint offsets
954  bands : `list` [`str`]
955  List of band names from FGCM output
956  datasetConfig : `lsst.meas.algorithms.DatasetConfig`
957  Config for reference dataset
958  """
959 
960  self.log.info("Outputting standard stars to %s" % (datasetConfig.ref_dataset_name))
961 
962  indexer = IndexerRegistry[self.config.datasetConfig.indexer.name](
963  self.config.datasetConfig.indexer.active)
964 
965  # We determine the conversion from the native units (typically radians) to
966  # degrees for the first star. This allows us to treat coord_ra/coord_dec as
967  # numpy arrays rather than Angles, which would we approximately 600x slower.
968  # TODO: Fix this after DM-16524 (HtmIndexer.indexPoints should take coords
969  # (as Angles) for input
970  conv = stdCat[0]['coord_ra'].asDegrees()/float(stdCat[0]['coord_ra'])
971  indices = np.array(indexer.indexPoints(stdCat['coord_ra']*conv,
972  stdCat['coord_dec']*conv))
973 
974  formattedCat = self._formatCatalog(stdCat, offsets, bands)
975 
976  # Write the master schema
977  dataId = indexer.makeDataId('master_schema',
978  datasetConfig.ref_dataset_name)
979  masterCat = afwTable.SimpleCatalog(formattedCat.schema)
980  addRefCatMetadata(masterCat)
981  butler.put(masterCat, 'ref_cat', dataId=dataId)
982 
983  # Break up the pixels using a histogram
984  h, rev = esutil.stat.histogram(indices, rev=True)
985  gd, = np.where(h > 0)
986  selected = np.zeros(len(formattedCat), dtype=bool)
987  for i in gd:
988  i1a = rev[rev[i]: rev[i + 1]]
989 
990  # the formattedCat afwTable can only be indexed with boolean arrays,
991  # and not numpy index arrays (see DM-16497). This little trick
992  # converts the index array into a boolean array
993  selected[:] = False
994  selected[i1a] = True
995 
996  # Write the individual pixel
997  dataId = indexer.makeDataId(indices[i1a[0]],
998  datasetConfig.ref_dataset_name)
999  butler.put(formattedCat[selected], 'ref_cat', dataId=dataId)
1000 
1001  # And save the dataset configuration
1002  dataId = indexer.makeDataId(None, datasetConfig.ref_dataset_name)
1003  butler.put(datasetConfig, 'ref_cat_config', dataId=dataId)
1004 
1005  self.log.info("Done outputting standard stars.")
1006 
1007  def _formatCatalog(self, fgcmStarCat, offsets, bands):
1008  """
1009  Turn an FGCM-formatted star catalog, applying zeropoint offsets.
1010 
1011  Parameters
1012  ----------
1013  fgcmStarCat : `lsst.afw.Table.SimpleCatalog`
1014  SimpleCatalog as output by fgcmcal
1015  offsets : `list` with len(self.bands) entries
1016  Zeropoint offsets to apply
1017  bands : `list` [`str`]
1018  List of band names from FGCM output
1019 
1020  Returns
1021  -------
1022  formattedCat: `lsst.afw.table.SimpleCatalog`
1023  SimpleCatalog suitable for using as a reference catalog
1024  """
1025 
1026  sourceMapper = afwTable.SchemaMapper(fgcmStarCat.schema)
1027  minSchema = LoadIndexedReferenceObjectsTask.makeMinimalSchema(bands,
1028  addCentroid=False,
1029  addIsResolved=True,
1030  coordErrDim=0)
1031  sourceMapper.addMinimalSchema(minSchema)
1032  for band in bands:
1033  sourceMapper.editOutputSchema().addField('%s_nGood' % (band), type=np.int32)
1034  sourceMapper.editOutputSchema().addField('%s_nTotal' % (band), type=np.int32)
1035  sourceMapper.editOutputSchema().addField('%s_nPsfCandidate' % (band), type=np.int32)
1036 
1037  formattedCat = afwTable.SimpleCatalog(sourceMapper.getOutputSchema())
1038  formattedCat.reserve(len(fgcmStarCat))
1039  formattedCat.extend(fgcmStarCat, mapper=sourceMapper)
1040 
1041  # Note that we don't have to set `resolved` because the default is False
1042 
1043  for b, band in enumerate(bands):
1044  mag = fgcmStarCat['mag_std_noabs'][:, b].astype(np.float64) + offsets[b]
1045  # We want fluxes in nJy from calibrated AB magnitudes
1046  # (after applying offset). Updated after RFC-549 and RFC-575.
1047  flux = (mag*units.ABmag).to_value(units.nJy)
1048  fluxErr = (np.log(10.)/2.5)*flux*fgcmStarCat['magErr_std'][:, b].astype(np.float64)
1049 
1050  formattedCat['%s_flux' % (band)][:] = flux
1051  formattedCat['%s_fluxErr' % (band)][:] = fluxErr
1052  formattedCat['%s_nGood' % (band)][:] = fgcmStarCat['ngood'][:, b]
1053  formattedCat['%s_nTotal' % (band)][:] = fgcmStarCat['ntotal'][:, b]
1054  formattedCat['%s_nPsfCandidate' % (band)][:] = fgcmStarCat['npsfcand'][:, b]
1055 
1056  addRefCatMetadata(formattedCat)
1057 
1058  return formattedCat
1059 
1060  def _outputZeropoints(self, camera, zptCat, visitCat, offsets, bands,
1061  physicalFilterMap, returnCatalogs=True,
1062  tract=None):
1063  """Output the zeropoints in fgcm_photoCalib format.
1064 
1065  Parameters
1066  ----------
1067  camera : `lsst.afw.cameraGeom.Camera`
1068  Camera from the butler.
1069  zptCat : `lsst.afw.table.BaseCatalog`
1070  FGCM zeropoint catalog from `FgcmFitCycleTask`.
1071  visitCat : `lsst.afw.table.BaseCatalog`
1072  FGCM visitCat from `FgcmBuildStarsTask`.
1073  offsets : `numpy.array`
1074  Float array of absolute calibration offsets, one for each filter.
1075  bands : `list` [`str`]
1076  List of band names from FGCM output.
1077  physicalFilterMap : `dict`
1078  Dictionary of mappings from physical filter to FGCM band.
1079  returnCatalogs : `bool`, optional
1080  Return photoCalibs as per-visit exposure catalogs.
1081  tract: `int`, optional
1082  Tract number to output. Default is None (global calibration)
1083 
1084  Returns
1085  -------
1086  photoCalibs : `generator` [(`int`, `int`, `str`, `lsst.afw.image.PhotoCalib`)]
1087  Generator that returns (visit, ccd, filtername, photoCalib) tuples.
1088  (returned if returnCatalogs is False).
1089  photoCalibCatalogs : `generator` [(`int`, `lsst.afw.table.ExposureCatalog`)]
1090  Generator that returns (visit, exposureCatalog) tuples.
1091  (returned if returnCatalogs is True).
1092  """
1093  # Select visit/ccds where we have a calibration
1094  # This includes ccds where we were able to interpolate from neighboring
1095  # ccds.
1096  cannot_compute = fgcm.fgcmUtilities.zpFlagDict['CANNOT_COMPUTE_ZEROPOINT']
1097  selected = (((zptCat['fgcmFlag'] & cannot_compute) == 0)
1098  & (zptCat['fgcmZptVar'] > 0.0)
1099  & (zptCat['fgcmZpt'] > FGCM_ILLEGAL_VALUE))
1100 
1101  # Log warnings for any visit which has no valid zeropoints
1102  badVisits = np.unique(zptCat['visit'][~selected])
1103  goodVisits = np.unique(zptCat['visit'][selected])
1104  allBadVisits = badVisits[~np.isin(badVisits, goodVisits)]
1105  for allBadVisit in allBadVisits:
1106  self.log.warning(f'No suitable photoCalib for visit {allBadVisit}')
1107 
1108  # Get a mapping from filtername to the offsets
1109  offsetMapping = {}
1110  for f in physicalFilterMap:
1111  # Not every filter in the map will necesarily have a band.
1112  if physicalFilterMap[f] in bands:
1113  offsetMapping[f] = offsets[bands.index(physicalFilterMap[f])]
1114 
1115  # Get a mapping from "ccd" to the ccd index used for the scaling
1116  ccdMapping = {}
1117  for ccdIndex, detector in enumerate(camera):
1118  ccdMapping[detector.getId()] = ccdIndex
1119 
1120  # And a mapping to get the flat-field scaling values
1121  scalingMapping = {}
1122  for rec in visitCat:
1123  scalingMapping[rec['visit']] = rec['scaling']
1124 
1125  if self.config.doComposeWcsJacobian:
1126  approxPixelAreaFields = computeApproxPixelAreaFields(camera)
1127 
1128  # The zptCat is sorted by visit, which is useful
1129  lastVisit = -1
1130  zptVisitCatalog = None
1131 
1132  metadata = dafBase.PropertyList()
1133  metadata.add("COMMENT", "Catalog id is detector id, sorted.")
1134  metadata.add("COMMENT", "Only detectors with data have entries.")
1135 
1136  for rec in zptCat[selected]:
1137  # Retrieve overall scaling
1138  scaling = scalingMapping[rec['visit']][ccdMapping[rec['detector']]]
1139 
1140  # The postCalibrationOffset describe any zeropoint offsets
1141  # to apply after the fgcm calibration. The first part comes
1142  # from the reference catalog match (used in testing). The
1143  # second part comes from the mean chromatic correction
1144  # (if configured).
1145  postCalibrationOffset = offsetMapping[rec['filtername']]
1146  if self.config.doApplyMeanChromaticCorrection:
1147  postCalibrationOffset += rec['fgcmDeltaChrom']
1148 
1149  fgcmSuperStarField = self._getChebyshevBoundedField(rec['fgcmfZptSstarCheb'],
1150  rec['fgcmfZptChebXyMax'])
1151  # Convert from FGCM AB to nJy
1152  fgcmZptField = self._getChebyshevBoundedField((rec['fgcmfZptCheb']*units.AB).to_value(units.nJy),
1153  rec['fgcmfZptChebXyMax'],
1154  offset=postCalibrationOffset,
1155  scaling=scaling)
1156 
1157  if self.config.doComposeWcsJacobian:
1158 
1159  fgcmField = afwMath.ProductBoundedField([approxPixelAreaFields[rec['detector']],
1160  fgcmSuperStarField,
1161  fgcmZptField])
1162  else:
1163  # The photoCalib is just the product of the fgcmSuperStarField and the
1164  # fgcmZptField
1165  fgcmField = afwMath.ProductBoundedField([fgcmSuperStarField, fgcmZptField])
1166 
1167  # The "mean" calibration will be set to the center of the ccd for reference
1168  calibCenter = fgcmField.evaluate(fgcmField.getBBox().getCenter())
1169  calibErr = (np.log(10.0)/2.5)*calibCenter*np.sqrt(rec['fgcmZptVar'])
1170  photoCalib = afwImage.PhotoCalib(calibrationMean=calibCenter,
1171  calibrationErr=calibErr,
1172  calibration=fgcmField,
1173  isConstant=False)
1174 
1175  if not returnCatalogs:
1176  # Return individual photoCalibs
1177  yield (int(rec['visit']), int(rec['detector']), rec['filtername'], photoCalib)
1178  else:
1179  # Return full per-visit exposure catalogs
1180  if rec['visit'] != lastVisit:
1181  # This is a new visit. If the last visit was not -1, yield
1182  # the ExposureCatalog
1183  if lastVisit > -1:
1184  # ensure that the detectors are in sorted order, for fast lookups
1185  zptVisitCatalog.sort()
1186  yield (int(lastVisit), zptVisitCatalog)
1187  else:
1188  # We need to create a new schema
1189  zptExpCatSchema = afwTable.ExposureTable.makeMinimalSchema()
1190  zptExpCatSchema.addField('visit', type='I', doc='Visit number')
1191 
1192  # And start a new one
1193  zptVisitCatalog = afwTable.ExposureCatalog(zptExpCatSchema)
1194  zptVisitCatalog.setMetadata(metadata)
1195 
1196  lastVisit = int(rec['visit'])
1197 
1198  catRecord = zptVisitCatalog.addNew()
1199  catRecord['id'] = int(rec['detector'])
1200  catRecord['visit'] = rec['visit']
1201  catRecord.setPhotoCalib(photoCalib)
1202 
1203  # Final output of last exposure catalog
1204  if returnCatalogs:
1205  # ensure that the detectors are in sorted order, for fast lookups
1206  zptVisitCatalog.sort()
1207  yield (int(lastVisit), zptVisitCatalog)
1208 
1209  def _getChebyshevBoundedField(self, coefficients, xyMax, offset=0.0, scaling=1.0):
1210  """
1211  Make a ChebyshevBoundedField from fgcm coefficients, with optional offset
1212  and scaling.
1213 
1214  Parameters
1215  ----------
1216  coefficients: `numpy.array`
1217  Flattened array of chebyshev coefficients
1218  xyMax: `list` of length 2
1219  Maximum x and y of the chebyshev bounding box
1220  offset: `float`, optional
1221  Absolute calibration offset. Default is 0.0
1222  scaling: `float`, optional
1223  Flat scaling value from fgcmBuildStars. Default is 1.0
1224 
1225  Returns
1226  -------
1227  boundedField: `lsst.afw.math.ChebyshevBoundedField`
1228  """
1229 
1230  orderPlus1 = int(np.sqrt(coefficients.size))
1231  pars = np.zeros((orderPlus1, orderPlus1))
1232 
1233  bbox = lsst.geom.Box2I(lsst.geom.Point2I(0.0, 0.0),
1234  lsst.geom.Point2I(*xyMax))
1235 
1236  pars[:, :] = (coefficients.reshape(orderPlus1, orderPlus1)
1237  * (10.**(offset/-2.5))*scaling)
1238 
1239  boundedField = afwMath.ChebyshevBoundedField(bbox, pars)
1240 
1241  return boundedField
1242 
1243  def _outputAtmospheres(self, dataRefDict, atmCat):
1244  """
1245  Output the atmospheres.
1246 
1247  Parameters
1248  ----------
1249  dataRefDict : `dict`
1250  All dataRefs are `lsst.daf.persistence.ButlerDataRef` (gen2) or
1251  `lsst.daf.butler.DeferredDatasetHandle` (gen3)
1252  dataRef dictionary with keys:
1253 
1254  ``"fgcmLookUpTable"``
1255  dataRef for the FGCM look-up table.
1256  atmCat : `lsst.afw.table.BaseCatalog`
1257  FGCM atmosphere parameter catalog from fgcmFitCycleTask.
1258 
1259  Returns
1260  -------
1261  atmospheres : `generator` [(`int`, `lsst.afw.image.TransmissionCurve`)]
1262  Generator that returns (visit, transmissionCurve) tuples.
1263  """
1264  # First, we need to grab the look-up table and key info
1265  lutCat = dataRefDict['fgcmLookUpTable'].get()
1266 
1267  atmosphereTableName = lutCat[0]['tablename']
1268  elevation = lutCat[0]['elevation']
1269  atmLambda = lutCat[0]['atmLambda']
1270  lutCat = None
1271 
1272  # Make the atmosphere table if possible
1273  try:
1274  atmTable = fgcm.FgcmAtmosphereTable.initWithTableName(atmosphereTableName)
1275  atmTable.loadTable()
1276  except IOError:
1277  atmTable = None
1278 
1279  if atmTable is None:
1280  # Try to use MODTRAN instead
1281  try:
1282  modGen = fgcm.ModtranGenerator(elevation)
1283  lambdaRange = np.array([atmLambda[0], atmLambda[-1]])/10.
1284  lambdaStep = (atmLambda[1] - atmLambda[0])/10.
1285  except (ValueError, IOError) as e:
1286  raise RuntimeError("FGCM look-up-table generated with modtran, "
1287  "but modtran not configured to run.") from e
1288 
1289  zenith = np.degrees(np.arccos(1./atmCat['secZenith']))
1290 
1291  for i, visit in enumerate(atmCat['visit']):
1292  if atmTable is not None:
1293  # Interpolate the atmosphere table
1294  atmVals = atmTable.interpolateAtmosphere(pmb=atmCat[i]['pmb'],
1295  pwv=atmCat[i]['pwv'],
1296  o3=atmCat[i]['o3'],
1297  tau=atmCat[i]['tau'],
1298  alpha=atmCat[i]['alpha'],
1299  zenith=zenith[i],
1300  ctranslamstd=[atmCat[i]['cTrans'],
1301  atmCat[i]['lamStd']])
1302  else:
1303  # Run modtran
1304  modAtm = modGen(pmb=atmCat[i]['pmb'],
1305  pwv=atmCat[i]['pwv'],
1306  o3=atmCat[i]['o3'],
1307  tau=atmCat[i]['tau'],
1308  alpha=atmCat[i]['alpha'],
1309  zenith=zenith[i],
1310  lambdaRange=lambdaRange,
1311  lambdaStep=lambdaStep,
1312  ctranslamstd=[atmCat[i]['cTrans'],
1313  atmCat[i]['lamStd']])
1314  atmVals = modAtm['COMBINED']
1315 
1316  # Now need to create something to persist...
1317  curve = TransmissionCurve.makeSpatiallyConstant(throughput=atmVals,
1318  wavelengths=atmLambda,
1319  throughputAtMin=atmVals[0],
1320  throughputAtMax=atmVals[-1])
1321 
1322  yield (int(visit), curve)
int min
A group of labels for a filter in an exposure or coadd.
Definition: FilterLabel.h:58
The photometric calibration of an exposure.
Definition: PhotoCalib.h:114
A BoundedField based on 2-d Chebyshev polynomials of the first kind.
A BoundedField that lazily multiplies a sequence of other BoundedFields.
Custom catalog class for ExposureRecord/Table.
Definition: Exposure.h:311
Defines the fields and offsets for a table.
Definition: Schema.h:51
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
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
An integer coordinate rectangle.
Definition: Box.h:55
daf::base::PropertyList * list
Definition: fits.cc:913
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
def computeApproxPixelAreaFields(camera)
Definition: utilities.py:499
def run(self, coaddExposures, bbox, wcs)
Definition: getTemplate.py:603
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.