23 """Make the final fgcmcal output products.
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
41 from astropy
import units
58 from lsst.utils.timer
import timeMethod
60 from .utilities
import computeApproxPixelAreaFields
61 from .utilities
import lookupStaticCalibrations
62 from .utilities
import FGCM_ILLEGAL_VALUE
66 __all__ = [
'FgcmOutputProductsConfig',
'FgcmOutputProductsTask',
'FgcmOutputProductsRunner']
70 dimensions=(
"instrument",),
71 defaultTemplates={
"cycleNumber":
"0"}):
72 camera = connectionTypes.PrerequisiteInput(
73 doc=
"Camera instrument",
75 storageClass=
"Camera",
76 dimensions=(
"instrument",),
77 lookupFunction=lookupStaticCalibrations,
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",),
90 fgcmVisitCatalog = connectionTypes.Input(
91 doc=
"Catalog of visit information for fgcm",
92 name=
"fgcmVisitCatalog",
93 storageClass=
"Catalog",
94 dimensions=(
"instrument",),
98 fgcmStandardStars = connectionTypes.Input(
99 doc=
"Catalog of standard star data from fgcm fit",
100 name=
"fgcmStandardStars{cycleNumber}",
101 storageClass=
"SimpleCatalog",
102 dimensions=(
"instrument",),
106 fgcmZeropoints = connectionTypes.Input(
107 doc=
"Catalog of zeropoints from fgcm fit",
108 name=
"fgcmZeropoints{cycleNumber}",
109 storageClass=
"Catalog",
110 dimensions=(
"instrument",),
114 fgcmAtmosphereParameters = connectionTypes.Input(
115 doc=
"Catalog of atmosphere parameters from fgcm fit",
116 name=
"fgcmAtmosphereParameters{cycleNumber}",
117 storageClass=
"Catalog",
118 dimensions=(
"instrument",),
122 refCat = connectionTypes.PrerequisiteInput(
123 doc=
"Reference catalog to use for photometric calibration",
125 storageClass=
"SimpleCatalog",
126 dimensions=(
"skypix",),
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",),
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",
150 fgcmOffsets = connectionTypes.Output(
151 doc=
"Per-band offsets computed from doReferenceCalibration",
152 name=
"fgcmReferenceCalibrationOffsets",
153 storageClass=
"Catalog",
154 dimensions=(
"instrument",),
158 def __init__(self, *, config=None):
159 super().__init__(config=config)
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")
166 if config.doRefcatOutput:
167 raise ValueError(
"FgcmOutputProductsTask (Gen3) does not support doRefcatOutput")
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")
179 class FgcmOutputProductsConfig(pipeBase.PipelineTaskConfig,
180 pipelineConnections=FgcmOutputProductsConnections):
181 """Config for FgcmOutputProductsTask"""
183 cycleNumber = pexConfig.Field(
184 doc=
"Final fit cycle from FGCM fit",
188 physicalFilterMap = pexConfig.DictField(
189 doc=
"Mapping from 'physicalFilter' to band.",
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."),
203 doRefcatOutput = pexConfig.Field(
204 doc=
"Output standard stars in reference catalog format",
208 doAtmosphereOutput = pexConfig.Field(
209 doc=
"Output atmospheres in transmission_atmosphere_fgcm format",
213 doZeropointOutput = pexConfig.Field(
214 doc=
"Output zeropoints in fgcm_photoCalib format",
218 doComposeWcsJacobian = pexConfig.Field(
219 doc=
"Compose Jacobian of WCS with fgcm calibration for output photoCalib?",
223 doApplyMeanChromaticCorrection = pexConfig.Field(
224 doc=
"Apply the mean chromatic correction to the zeropoints?",
228 refObjLoader = pexConfig.ConfigurableField(
229 target=LoadIndexedReferenceObjectsTask,
230 doc=
"reference object loader for 'absolute' photometric calibration",
232 photoCal = pexConfig.ConfigurableField(
234 doc=
"task to perform 'absolute' calibration",
236 referencePixelizationNside = pexConfig.Field(
237 doc=
"Healpix nside to pixelize catalog to compare to reference catalog",
241 referencePixelizationMinStars = pexConfig.Field(
242 doc=(
"Minimum number of stars per healpix pixel to select for comparison"
243 "to the specified reference catalog"),
247 referenceMinMatch = pexConfig.Field(
248 doc=
"Minimum number of stars matched to reference catalog to be used in statistics",
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."),
260 datasetConfig = pexConfig.ConfigField(
262 doc=
"Configuration for writing/reading ingested catalog",
266 pexConfig.Config.setDefaults(self)
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
296 self.connections.cycleNumber = str(self.cycleNumber)
299 class FgcmOutputProductsRunner(pipeBase.ButlerInitializedTaskRunner):
300 """Subclass of TaskRunner for fgcmOutputProductsTask
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.
308 def getTargetList(parsedCmd):
310 Return a list with one element, the butler.
312 return [parsedCmd.butler]
314 def __call__(self, butler):
318 butler: `lsst.daf.persistence.Butler`
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)
327 task = self.TaskClass(butler=butler, config=self.config, log=self.log)
331 results = task.runDataRef(butler)
334 results = task.runDataRef(butler)
335 except Exception
as e:
337 task.log.fatal(
"Failed: %s" % e)
338 if not isinstance(e, pipeBase.TaskError):
339 traceback.print_exc(file=sys.stderr)
341 task.writeMetadata(butler)
343 if self.doReturnResults:
345 return [pipeBase.Struct(exitStatus=exitStatus,
348 return [pipeBase.Struct(exitStatus=exitStatus)]
350 def run(self, parsedCmd):
352 Run the task, with no multiprocessing
356 parsedCmd: `lsst.pipe.base.ArgumentParser` parsed command line
361 if self.precall(parsedCmd):
362 targetList = self.getTargetList(parsedCmd)
364 resultList = self(targetList[0])
369 class FgcmOutputProductsTask(pipeBase.PipelineTask, pipeBase.CmdLineTask):
371 Output products from FGCM global calibration.
374 ConfigClass = FgcmOutputProductsConfig
375 RunnerClass = FgcmOutputProductsRunner
376 _DefaultName =
"fgcmOutputProducts"
378 def __init__(self, butler=None, **kwargs):
379 super().__init__(**kwargs)
382 def _getMetadataName(self):
385 def runQuantum(self, butlerQC, inputRefs, outputRefs):
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)
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}
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}
402 if self.config.doReferenceCalibration:
403 refConfig = self.config.refObjLoader
405 for ref
in inputRefs.refCat],
406 refCats=butlerQC.get(inputRefs.refCat),
410 self.refObjLoader =
None
412 struct = self.run(dataRefDict, self.config.physicalFilterMap, returnCatalogs=
True)
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.")
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.")
428 if self.config.doReferenceCalibration:
431 schema.addField(
'offset', type=np.float64,
432 doc=
"Post-process calibration offset (mag)")
434 offsetCat.resize(len(struct.offsets))
435 offsetCat[
'offset'][:] = struct.offsets
437 butlerQC.put(offsetCat, outputRefs.fgcmOffsets)
442 def runDataRef(self, butler):
444 Make FGCM output products for use in the stack
448 butler: `lsst.daf.persistence.Butler`
450 Final fit cycle number, override config.
454 offsets: `lsst.pipe.base.Struct`
455 A structure with array of zeropoint offsets
460 Raised if any one of the following is true:
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`.
475 if self.config.doReferenceCalibration:
477 self.makeSubtask(
"refObjLoader", butler=butler)
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")
486 if butler.datasetExists(
'fgcmBuildStarsTable_config'):
487 fgcmBuildStarsConfig = butler.get(
'fgcmBuildStarsTable_config')
489 fgcmBuildStarsConfig = butler.get(
'fgcmBuildStars_config')
490 visitDataRefName = fgcmBuildStarsConfig.visitDataRefName
491 ccdDataRefName = fgcmBuildStarsConfig.ccdDataRefName
492 physicalFilterMap = fgcmBuildStarsConfig.physicalFilterMap
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.")
498 if not self.config.doComposeWcsJacobian
and fgcmBuildStarsConfig.doApplyWcsJacobian:
499 self.log.
warning(
"Jacobian was applied in build-stars but doComposeWcsJacobian is not set.")
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}.")
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))
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))
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)
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)
531 struct = self.run(dataRefDict, physicalFilterMap, butler=butler, returnCatalogs=
False)
533 if struct.photoCalibs
is not None:
534 self.log.
info(
"Outputting photoCalib files.")
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})
542 self.log.
info(
"Done outputting photoCalib files.")
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.")
551 return pipeBase.Struct(offsets=struct.offsets)
553 def run(self, dataRefDict, physicalFilterMap, returnCatalogs=True, butler=None):
554 """Run the output products task.
559 All dataRefs are `lsst.daf.persistence.ButlerDataRef` (gen2) or
560 `lsst.daf.butler.DeferredDatasetHandle` (gen3)
561 dataRef dictionary with keys:
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.
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
586 retStruct : `lsst.pipe.base.Struct`
587 Output structure with keys:
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).
600 stdCat = dataRefDict[
'fgcmStandardStars'].get()
601 md = stdCat.getMetadata()
602 bands = md.getArray(
'BANDS')
604 if self.config.doReferenceCalibration:
605 lutCat = dataRefDict[
'fgcmLookUpTable'].get()
606 offsets = self._computeReferenceOffsets(stdCat, lutCat, physicalFilterMap, bands)
608 offsets = np.zeros(len(bands))
611 if self.config.doRefcatOutput
and butler
is not None:
612 self._outputStandardStars(butler, stdCat, offsets, bands, self.config.datasetConfig)
616 if self.config.doZeropointOutput:
617 zptCat = dataRefDict[
'fgcmZeropoints'].get()
618 visitCat = dataRefDict[
'fgcmVisitCatalog'].get()
620 pcgen = self._outputZeropoints(dataRefDict[
'camera'], zptCat, visitCat, offsets, bands,
621 physicalFilterMap, returnCatalogs=returnCatalogs)
625 if self.config.doAtmosphereOutput:
626 atmCat = dataRefDict[
'fgcmAtmosphereParameters'].get()
627 atmgen = self._outputAtmospheres(dataRefDict, atmCat)
631 retStruct = pipeBase.Struct(offsets=offsets,
634 retStruct.photoCalibCatalogs = pcgen
636 retStruct.photoCalibs = pcgen
640 def generateTractOutputProducts(self, dataRefDict, tract,
641 visitCat, zptCat, atmCat, stdCat,
642 fgcmBuildStarsConfig,
646 Generate the output products for a given tract, as specified in the config.
648 This method is here to have an alternate entry-point for
654 All dataRefs are `lsst.daf.persistence.ButlerDataRef` (gen2) or
655 `lsst.daf.butler.DeferredDatasetHandle` (gen3)
656 dataRef dictionary with keys:
659 Camera object (`lsst.afw.cameraGeom.Camera`)
660 ``"fgcmLookUpTable"``
661 dataRef for the FGCM look-up table.
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
681 retStruct : `lsst.pipe.base.Struct`
682 Output structure with keys:
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).
695 physicalFilterMap = fgcmBuildStarsConfig.physicalFilterMap
697 md = stdCat.getMetadata()
698 bands = md.getArray(
'BANDS')
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.")
704 if not self.config.doComposeWcsJacobian
and fgcmBuildStarsConfig.doApplyWcsJacobian:
705 self.log.
warning(
"Jacobian was applied in build-stars but doComposeWcsJacobian is not set.")
707 if self.config.doReferenceCalibration:
708 lutCat = dataRefDict[
'fgcmLookUpTable'].get()
709 offsets = self._computeReferenceOffsets(stdCat, lutCat, bands, physicalFilterMap)
711 offsets = np.zeros(len(bands))
713 if self.config.doRefcatOutput
and butler
is not None:
715 datasetConfig = copy.copy(self.config.datasetConfig)
716 datasetConfig.ref_dataset_name =
'%s_%d' % (self.config.datasetConfig.ref_dataset_name,
718 self._outputStandardStars(butler, stdCat, offsets, bands, datasetConfig)
720 if self.config.doZeropointOutput:
721 pcgen = self._outputZeropoints(dataRefDict[
'camera'], zptCat, visitCat, offsets, bands,
722 physicalFilterMap, returnCatalogs=returnCatalogs)
726 if self.config.doAtmosphereOutput:
727 atmgen = self._outputAtmospheres(dataRefDict, atmCat)
731 retStruct = pipeBase.Struct(offsets=offsets,
734 retStruct.photoCalibCatalogs = pcgen
736 retStruct.photoCalibs = pcgen
740 def _computeReferenceOffsets(self, stdCat, lutCat, physicalFilterMap, bands):
742 Compute offsets relative to a reference catalog.
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
752 stdCat : `lsst.afw.table.SimpleCatalog`
754 lutCat : `lsst.afw.table.SimpleCatalog`
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
762 offsets : `numpy.array` of floats
763 Per band zeropoint offsets
769 minObs = stdCat[
'ngood'].
min(axis=1)
771 goodStars = (minObs >= 1)
772 stdCat = stdCat[goodStars]
774 self.log.
info(
"Found %d stars with at least 1 good observation in each band" %
781 lutPhysicalFilters = lutCat[0][
'physicalFilters'].split(
',')
782 lutStdPhysicalFilters = lutCat[0][
'stdPhysicalFilters'].split(
',')
783 physicalFilterMapBands =
list(physicalFilterMap.values())
784 physicalFilterMapFilters =
list(physicalFilterMap.keys())
788 physicalFilterMapIndex = physicalFilterMapBands.index(band)
789 physicalFilter = physicalFilterMapFilters[physicalFilterMapIndex]
791 lutPhysicalFilterIndex = lutPhysicalFilters.index(physicalFilter)
792 stdPhysicalFilter = lutStdPhysicalFilters[lutPhysicalFilterIndex]
794 physical=stdPhysicalFilter))
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',
817 theta = np.pi/2. - stdCat[
'coord_dec']
818 phi = stdCat[
'coord_ra']
820 ipring = hp.ang2pix(self.config.referencePixelizationNside, theta, phi)
821 h, rev = esutil.stat.histogram(ipring, rev=
True)
823 gdpix, = np.where(h >= self.config.referencePixelizationMinStars)
825 self.log.
info(
"Found %d pixels (nside=%d) with at least %d good stars" %
827 self.config.referencePixelizationNside,
828 self.config.referencePixelizationMinStars))
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))
835 gdpix = np.random.choice(gdpix, size=self.config.referencePixelizationNPixels, replace=
False)
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]]]
845 selected = np.zeros(len(stdCat), dtype=bool)
847 refFluxFields = [
None]*len(bands)
849 for p_index, pix
in enumerate(gdpix):
850 i1a = rev[rev[pix]: rev[pix + 1]]
858 for b_index, filterLabel
in enumerate(filterLabels):
859 struct = self._computeOffsetOneBand(sourceMapper, badStarKey, b_index,
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
868 offsets = np.zeros(len(bands))
870 for b_index, band
in enumerate(bands):
872 ok, = np.where(results[
'nmatch'][:, b_index] >= self.config.referenceMinMatch)
873 offsets[b_index] = np.median(results[
'zp'][ok, b_index])
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)
882 def _computeOffsetOneBand(self, sourceMapper, badStarKey,
883 b_index, filterLabel, stdCat, selected, refFluxFields):
885 Compute the zeropoint offset between the fgcm stdCat and the reference
886 stars for one pixel in one band
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
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`
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
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'])
915 badStar = (stdCat[
'mag_std_noabs'][selected, b_index] > 90.0)
916 for rec
in sourceCat[badStar]:
917 rec.set(badStarKey,
True)
919 exposure = afwImage.ExposureF()
920 exposure.setFilterLabel(filterLabel)
922 if refFluxFields[b_index]
is None:
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
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,
936 schema=sourceCat.getSchema())
938 struct = calTask.run(exposure, sourceCat)
942 def _outputStandardStars(self, butler, stdCat, offsets, bands, datasetConfig):
944 Output standard stars in indexed reference catalog format.
945 This is not currently supported in Gen3.
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
960 self.log.
info(
"Outputting standard stars to %s" % (datasetConfig.ref_dataset_name))
962 indexer = IndexerRegistry[self.config.datasetConfig.indexer.name](
963 self.config.datasetConfig.indexer.active)
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))
974 formattedCat = self._formatCatalog(stdCat, offsets, bands)
977 dataId = indexer.makeDataId(
'master_schema',
978 datasetConfig.ref_dataset_name)
981 butler.put(masterCat,
'ref_cat', dataId=dataId)
984 h, rev = esutil.stat.histogram(indices, rev=
True)
985 gd, = np.where(h > 0)
986 selected = np.zeros(len(formattedCat), dtype=bool)
988 i1a = rev[rev[i]: rev[i + 1]]
997 dataId = indexer.makeDataId(indices[i1a[0]],
998 datasetConfig.ref_dataset_name)
999 butler.put(formattedCat[selected],
'ref_cat', dataId=dataId)
1002 dataId = indexer.makeDataId(
None, datasetConfig.ref_dataset_name)
1003 butler.put(datasetConfig,
'ref_cat_config', dataId=dataId)
1005 self.log.
info(
"Done outputting standard stars.")
1007 def _formatCatalog(self, fgcmStarCat, offsets, bands):
1009 Turn an FGCM-formatted star catalog, applying zeropoint offsets.
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
1022 formattedCat: `lsst.afw.table.SimpleCatalog`
1023 SimpleCatalog suitable for using as a reference catalog
1027 minSchema = LoadIndexedReferenceObjectsTask.makeMinimalSchema(bands,
1031 sourceMapper.addMinimalSchema(minSchema)
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)
1038 formattedCat.reserve(len(fgcmStarCat))
1039 formattedCat.extend(fgcmStarCat, mapper=sourceMapper)
1043 for b, band
in enumerate(bands):
1044 mag = fgcmStarCat[
'mag_std_noabs'][:, b].astype(np.float64) + offsets[b]
1047 flux = (mag*units.ABmag).to_value(units.nJy)
1048 fluxErr = (np.log(10.)/2.5)*flux*fgcmStarCat[
'magErr_std'][:, b].astype(np.float64)
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]
1060 def _outputZeropoints(self, camera, zptCat, visitCat, offsets, bands,
1061 physicalFilterMap, returnCatalogs=True,
1063 """Output the zeropoints in fgcm_photoCalib format.
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)
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).
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))
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}')
1110 for f
in physicalFilterMap:
1112 if physicalFilterMap[f]
in bands:
1113 offsetMapping[f] = offsets[bands.index(physicalFilterMap[f])]
1117 for ccdIndex, detector
in enumerate(camera):
1118 ccdMapping[detector.getId()] = ccdIndex
1122 for rec
in visitCat:
1123 scalingMapping[rec[
'visit']] = rec[
'scaling']
1125 if self.config.doComposeWcsJacobian:
1130 zptVisitCatalog =
None
1133 metadata.add(
"COMMENT",
"Catalog id is detector id, sorted.")
1134 metadata.add(
"COMMENT",
"Only detectors with data have entries.")
1136 for rec
in zptCat[selected]:
1138 scaling = scalingMapping[rec[
'visit']][ccdMapping[rec[
'detector']]]
1145 postCalibrationOffset = offsetMapping[rec[
'filtername']]
1146 if self.config.doApplyMeanChromaticCorrection:
1147 postCalibrationOffset += rec[
'fgcmDeltaChrom']
1149 fgcmSuperStarField = self._getChebyshevBoundedField(rec[
'fgcmfZptSstarCheb'],
1150 rec[
'fgcmfZptChebXyMax'])
1152 fgcmZptField = self._getChebyshevBoundedField((rec[
'fgcmfZptCheb']*units.AB).to_value(units.nJy),
1153 rec[
'fgcmfZptChebXyMax'],
1154 offset=postCalibrationOffset,
1157 if self.config.doComposeWcsJacobian:
1168 calibCenter = fgcmField.evaluate(fgcmField.getBBox().getCenter())
1169 calibErr = (np.log(10.0)/2.5)*calibCenter*np.sqrt(rec[
'fgcmZptVar'])
1171 calibrationErr=calibErr,
1172 calibration=fgcmField,
1175 if not returnCatalogs:
1177 yield (int(rec[
'visit']), int(rec[
'detector']), rec[
'filtername'], photoCalib)
1180 if rec[
'visit'] != lastVisit:
1185 zptVisitCatalog.sort()
1186 yield (int(lastVisit), zptVisitCatalog)
1189 zptExpCatSchema = afwTable.ExposureTable.makeMinimalSchema()
1190 zptExpCatSchema.addField(
'visit', type=
'I', doc=
'Visit number')
1194 zptVisitCatalog.setMetadata(metadata)
1196 lastVisit = int(rec[
'visit'])
1198 catRecord = zptVisitCatalog.addNew()
1199 catRecord[
'id'] = int(rec[
'detector'])
1200 catRecord[
'visit'] = rec[
'visit']
1201 catRecord.setPhotoCalib(photoCalib)
1206 zptVisitCatalog.sort()
1207 yield (int(lastVisit), zptVisitCatalog)
1209 def _getChebyshevBoundedField(self, coefficients, xyMax, offset=0.0, scaling=1.0):
1211 Make a ChebyshevBoundedField from fgcm coefficients, with optional offset
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
1227 boundedField: `lsst.afw.math.ChebyshevBoundedField`
1230 orderPlus1 = int(np.sqrt(coefficients.size))
1231 pars = np.zeros((orderPlus1, orderPlus1))
1236 pars[:, :] = (coefficients.reshape(orderPlus1, orderPlus1)
1237 * (10.**(offset/-2.5))*scaling)
1243 def _outputAtmospheres(self, dataRefDict, atmCat):
1245 Output the atmospheres.
1249 dataRefDict : `dict`
1250 All dataRefs are `lsst.daf.persistence.ButlerDataRef` (gen2) or
1251 `lsst.daf.butler.DeferredDatasetHandle` (gen3)
1252 dataRef dictionary with keys:
1254 ``"fgcmLookUpTable"``
1255 dataRef for the FGCM look-up table.
1256 atmCat : `lsst.afw.table.BaseCatalog`
1257 FGCM atmosphere parameter catalog from fgcmFitCycleTask.
1261 atmospheres : `generator` [(`int`, `lsst.afw.image.TransmissionCurve`)]
1262 Generator that returns (visit, transmissionCurve) tuples.
1265 lutCat = dataRefDict[
'fgcmLookUpTable'].get()
1267 atmosphereTableName = lutCat[0][
'tablename']
1268 elevation = lutCat[0][
'elevation']
1269 atmLambda = lutCat[0][
'atmLambda']
1274 atmTable = fgcm.FgcmAtmosphereTable.initWithTableName(atmosphereTableName)
1275 atmTable.loadTable()
1279 if atmTable
is None:
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
1289 zenith = np.degrees(np.arccos(1./atmCat[
'secZenith']))
1291 for i, visit
in enumerate(atmCat[
'visit']):
1292 if atmTable
is not None:
1294 atmVals = atmTable.interpolateAtmosphere(pmb=atmCat[i][
'pmb'],
1295 pwv=atmCat[i][
'pwv'],
1297 tau=atmCat[i][
'tau'],
1298 alpha=atmCat[i][
'alpha'],
1300 ctranslamstd=[atmCat[i][
'cTrans'],
1301 atmCat[i][
'lamStd']])
1304 modAtm = modGen(pmb=atmCat[i][
'pmb'],
1305 pwv=atmCat[i][
'pwv'],
1307 tau=atmCat[i][
'tau'],
1308 alpha=atmCat[i][
'alpha'],
1310 lambdaRange=lambdaRange,
1311 lambdaStep=lambdaStep,
1312 ctranslamstd=[atmCat[i][
'cTrans'],
1313 atmCat[i][
'lamStd']])
1314 atmVals = modAtm[
'COMBINED']
1317 curve = TransmissionCurve.makeSpatiallyConstant(throughput=atmVals,
1318 wavelengths=atmLambda,
1319 throughputAtMin=atmVals[0],
1320 throughputAtMax=atmVals[-1])
1322 yield (int(visit), curve)
A group of labels for a filter in an exposure or coadd.
The photometric calibration of an exposure.
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.
Defines the fields and offsets for a table.
A mapping between the keys of two Schemas, used to copy data between them.
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
Class for storing ordered metadata with comments.
An integer coordinate rectangle.
daf::base::PropertyList * list
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
def computeApproxPixelAreaFields(camera)
def run(self, coaddExposures, bbox, wcs)
def addRefCatMetadata(catalog)
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.