LSST Applications g013ef56533+7c9321ec0f,g042eb84c57+c6cfa41bc3,g199a45376c+0ba108daf9,g1fd858c14a+fcad0d0313,g210f2d0738+c0f94c6586,g262e1987ae+a7e710680e,g29ae962dfc+fb55f2edb0,g2ac17093b6+61d6563b1e,g2b1d02342f+df6f932764,g2cef7863aa+aef1011c0b,g2f7ad74990+c0f94c6586,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53cf87ae69,g47891489e3+4316d04fff,g511e8cfd20+baa56acf6c,g53246c7159+8c5ae1fdc5,g54cd7ddccb+fd7ad03fde,g64539dfbff+c0f94c6586,g67b6fd64d1+4316d04fff,g67fd3c3899+c0f94c6586,g6985122a63+4316d04fff,g74acd417e5+ca833bee28,g786e29fd12+668abc6043,g81db2e9a8d+b2ec8e584f,g87389fa792+8856018cbb,g89139ef638+4316d04fff,g8d7436a09f+0a24083b20,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+11eb8fd5b8,gbf99507273+8c5ae1fdc5,gcdda8b9158+e4c84c9d5c,gce8aa8abaa+8c5ae1fdc5,gd7ef33dd92+4316d04fff,gdab6d2f7ff+ca833bee28,ge410e46f29+4316d04fff,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.40
LSST Data Management Base Package
Loading...
Searching...
No Matches
fgcmMakeLut.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 a look-up-table (LUT) for FGCM calibration.
24
25This task computes a look-up-table for the range in expected atmosphere
26variation and variation in instrumental throughput (as tracked by the
27transmission_filter products). By pre-computing linearized integrals,
28the FGCM fit is orders of magnitude faster for stars with a broad range
29of colors and observing bands, yielding precision at the 1-2 mmag level.
30
31Computing a LUT requires running MODTRAN or with a pre-generated
32atmosphere table packaged with fgcm.
33"""
34
35import numpy as np
36
37import lsst.pex.config as pexConfig
38import lsst.pipe.base as pipeBase
39from lsst.pipe.base import connectionTypes
40import lsst.afw.table as afwTable
41import lsst.afw.cameraGeom as afwCameraGeom
42from lsst.afw.image import TransmissionCurve
43from .utilities import lookupStaticCalibrations, countDetectors
44from astropy.table import Table
45import astropy.units as units
46
47import fgcm
48
49__all__ = ['FgcmMakeLutParametersConfig', 'FgcmMakeLutConfig', 'FgcmMakeLutTask', 'SensorCorrectionTerms']
50
51
52class SensorCorrectionTerms(pexConfig.Config):
53 refLambda = pexConfig.Field(
54 doc="Reference wavelength for first-order correction terms.",
55 dtype=float,
56 optional=False,
57 )
58 correctionTermDict = pexConfig.DictField(
59 doc="Mapping of detector number to first-order correction term.",
60 keytype=int,
61 itemtype=float,
62 default={},
63 )
64
65
66class FgcmMakeLutConnections(pipeBase.PipelineTaskConnections,
67 dimensions=('instrument',),
68 defaultTemplates={}):
69 camera = connectionTypes.PrerequisiteInput(
70 doc="Camera instrument",
71 name="camera",
72 storageClass="Camera",
73 dimensions=("instrument",),
74 isCalibration=True,
75 )
76
77 transmission_optics = connectionTypes.PrerequisiteInput(
78 doc="Optics transmission curve information",
79 name="transmission_optics",
80 storageClass="TransmissionCurve",
81 dimensions=("instrument",),
82 isCalibration=True,
83 deferLoad=True,
84 )
85
86 transmission_sensor = connectionTypes.PrerequisiteInput(
87 doc="Sensor transmission curve information",
88 name="transmission_sensor",
89 storageClass="TransmissionCurve",
90 dimensions=("instrument", "detector",),
91 lookupFunction=lookupStaticCalibrations,
92 isCalibration=True,
93 deferLoad=True,
94 multiple=True,
95 )
96
97 transmission_filter = connectionTypes.PrerequisiteInput(
98 doc="Filter transmission curve information",
99 name="transmission_filter",
100 storageClass="TransmissionCurve",
101 dimensions=("band", "instrument", "physical_filter",),
102 lookupFunction=lookupStaticCalibrations,
103 isCalibration=True,
104 deferLoad=True,
105 multiple=True,
106 )
107
108 transmission_filter_detector = connectionTypes.PrerequisiteInput(
109 doc="Filter transmission curve per detector",
110 name="transmission_filter_detector",
111 storageClass="TransmissionCurve",
112 dimensions=("instrument", "physical_filter", "detector",),
113 lookupFunction=lookupStaticCalibrations,
114 isCalibration=True,
115 deferLoad=True,
116 multiple=True,
117 )
118
119 fgcmLookUpTable = connectionTypes.Output(
120 doc=("Atmosphere + instrument look-up-table for FGCM throughput and "
121 "chromatic corrections."),
122 name="fgcmLookUpTable",
123 storageClass="Catalog",
124 dimensions=("instrument",),
125 )
126
127 fgcmStandardAtmosphere = connectionTypes.Output(
128 doc="Standard atmosphere used for FGCM calibration.",
129 name="fgcm_standard_atmosphere",
130 storageClass="TransmissionCurve",
131 dimensions=("instrument",),
132 )
133
134 fgcmStandardPassbands = connectionTypes.Output(
135 doc="Standard passbands used for FGCM calibration.",
136 name="fgcm_standard_passband",
137 storageClass="TransmissionCurve",
138 dimensions=("instrument", "physical_filter"),
139 multiple=True,
140 )
141
142 standardPassbands = connectionTypes.Output(
143 doc="Standard passbands, in astropy table format.",
144 name="standard_passband",
145 storageClass="ArrowAstropy",
146 dimensions=("instrument", "band"),
147 multiple=True,
148 )
149
150 def __init__(self, *, config=None):
151 if not config.doOpticsTransmission:
152 del self.transmission_optics
153 if not config.doSensorTransmission:
154 del self.transmission_sensor
155 if not config.doFilterDetectorTransmission:
157 else:
158 del self.transmission_filter
159
160
161class FgcmMakeLutParametersConfig(pexConfig.Config):
162 """Config for parameters if atmosphereTableName not available"""
163 # TODO: When DM-16511 is done, it will be possible to get the
164 # telescope elevation directly from the camera.
165 elevation = pexConfig.Field(
166 doc="Telescope elevation (m)",
167 dtype=float,
168 default=None,
169 )
170 pmbRange = pexConfig.ListField(
171 doc=("Barometric Pressure range (millibar) "
172 "Recommended range depends on the site."),
173 dtype=float,
174 default=None,
175 )
176 pmbSteps = pexConfig.Field(
177 doc="Barometric Pressure number of steps",
178 dtype=int,
179 default=5,
180 )
181 pwvRange = pexConfig.ListField(
182 doc=("Precipitable Water Vapor range (mm) "
183 "Recommended range depends on the site."),
184 dtype=float,
185 default=None,
186 )
187 pwvSteps = pexConfig.Field(
188 doc="Precipitable Water Vapor number of steps",
189 dtype=int,
190 default=15,
191 )
192 o3Range = pexConfig.ListField(
193 doc="Ozone range (dob)",
194 dtype=float,
195 default=[220.0, 310.0],
196 )
197 o3Steps = pexConfig.Field(
198 doc="Ozone number of steps",
199 dtype=int,
200 default=3,
201 )
202 tauRange = pexConfig.ListField(
203 doc="Aerosol Optical Depth range (unitless)",
204 dtype=float,
205 default=[0.002, 0.35],
206 )
207 tauSteps = pexConfig.Field(
208 doc="Aerosol Optical Depth number of steps",
209 dtype=int,
210 default=11,
211 )
212 alphaRange = pexConfig.ListField(
213 doc="Aerosol alpha range (unitless)",
214 dtype=float,
215 default=[0.0, 2.0],
216 )
217 alphaSteps = pexConfig.Field(
218 doc="Aerosol alpha number of steps",
219 dtype=int,
220 default=9,
221 )
222 zenithRange = pexConfig.ListField(
223 doc="Zenith angle range (degree)",
224 dtype=float,
225 default=[0.0, 70.0],
226 )
227 zenithSteps = pexConfig.Field(
228 doc="Zenith angle number of steps",
229 dtype=int,
230 default=21,
231 )
232 # Note that the standard atmosphere parameters depend on the observatory
233 # and elevation, and so these should be set on a per-camera basis.
234 pmbStd = pexConfig.Field(
235 doc=("Standard Atmosphere pressure (millibar); "
236 "Recommended default depends on the site."),
237 dtype=float,
238 default=None,
239 )
240 pwvStd = pexConfig.Field(
241 doc=("Standard Atmosphere PWV (mm); "
242 "Recommended default depends on the site."),
243 dtype=float,
244 default=None,
245 )
246 o3Std = pexConfig.Field(
247 doc="Standard Atmosphere O3 (dob)",
248 dtype=float,
249 default=263.0,
250 )
251 tauStd = pexConfig.Field(
252 doc="Standard Atmosphere aerosol optical depth",
253 dtype=float,
254 default=0.03,
255 )
256 alphaStd = pexConfig.Field(
257 doc="Standard Atmosphere aerosol alpha",
258 dtype=float,
259 default=1.0,
260 )
261 airmassStd = pexConfig.Field(
262 doc=("Standard Atmosphere airmass; "
263 "Recommended default depends on the survey strategy."),
264 dtype=float,
265 default=None,
266 )
267 lambdaNorm = pexConfig.Field(
268 doc="Aerosol Optical Depth normalization wavelength (Angstrom)",
269 dtype=float,
270 default=7750.0,
271 )
272 lambdaStep = pexConfig.Field(
273 doc="Wavelength step for generating atmospheres (nm)",
274 dtype=float,
275 default=0.5,
276 )
277 lambdaRange = pexConfig.ListField(
278 doc="Wavelength range for LUT (Angstrom)",
279 dtype=float,
280 default=[3000.0, 11000.0],
281 )
282
283
284class FgcmMakeLutConfig(pipeBase.PipelineTaskConfig,
285 pipelineConnections=FgcmMakeLutConnections):
286 """Config for FgcmMakeLutTask"""
287 physicalFilters = pexConfig.ListField(
288 doc="List of physicalFilter labels to generate look-up table.",
289 dtype=str,
290 default=[],
291 )
292 stdPhysicalFilterOverrideMap = pexConfig.DictField(
293 doc=("Override mapping from physical filter labels to 'standard' physical "
294 "filter labels. The 'standard' physical filter defines the transmission "
295 "curve that the FGCM standard bandpass will be based on. "
296 "Any filter not listed here will be mapped to "
297 "itself (e.g. g->g or HSC-G->HSC-G). Use this override for cross-"
298 "filter calibration such as HSC-R->HSC-R2 and HSC-I->HSC-I2."),
299 keytype=str,
300 itemtype=str,
301 default={},
302 )
303 atmosphereTableName = pexConfig.Field(
304 doc="FGCM name or filename of precomputed atmospheres",
305 dtype=str,
306 default=None,
307 optional=True,
308 )
309 useScienceDetectors = pexConfig.Field(
310 doc="Only use science detectors in LUT?",
311 dtype=bool,
312 default=True,
313 )
314 doOpticsTransmission = pexConfig.Field(
315 doc="Include optics transmission?",
316 dtype=bool,
317 default=True,
318 )
319 doSensorTransmission = pexConfig.Field(
320 doc="Include sensor transmission?",
321 dtype=bool,
322 default=True,
323 )
324 doFilterDetectorTransmission = pexConfig.Field(
325 doc="Use filter transmissions that are specified per-detector, rather "
326 "than a constant or radially-dependent filter transmission?",
327 dtype=bool,
328 default=False,
329 )
330 sensorCorrectionTermDict = pexConfig.ConfigDictField(
331 doc="Mapping of filter name to sensor correction terms.",
332 keytype=str,
333 itemtype=SensorCorrectionTerms,
334 default={},
335 )
336 parameters = pexConfig.ConfigField(
337 doc="Atmosphere parameters (required if no atmosphereTableName)",
338 dtype=FgcmMakeLutParametersConfig,
339 default=None,
340 check=None)
341
342 def validate(self):
343 """
344 Validate the config parameters.
345
346 This method behaves differently from the parent validate in the case
347 that atmosphereTableName is set. In this case, the config values
348 for standard values, step sizes, and ranges are loaded
349 directly from the specified atmosphereTableName.
350 """
351 # check that filterNames and stdFilterNames are okay
352 self._fields['physicalFilters'].validate(self)
353 self._fields['stdPhysicalFilterOverrideMap'].validate(self)
354
355 if self.atmosphereTableName is None:
356 # Validate the parameters
357 self._fields['parameters'].validate(self)
358
359
360class FgcmMakeLutTask(pipeBase.PipelineTask):
361 """
362 Make Look-Up Table for FGCM.
363
364 This task computes a look-up-table for the range in expected atmosphere
365 variation and variation in instrumental throughput (as tracked by the
366 transmission_filter products). By pre-computing linearized integrals,
367 the FGCM fit is orders of magnitude faster for stars with a broad range
368 of colors and observing bands, yielding precision at the 1-2 mmag level.
369
370 Computing a LUT requires running MODTRAN or with a pre-generated
371 atmosphere table packaged with fgcm.
372 """
373
374 ConfigClass = FgcmMakeLutConfig
375 _DefaultName = "fgcmMakeLut"
376
377 def __init__(self, initInputs=None, **kwargs):
378 super().__init__(**kwargs)
379
380 def runQuantum(self, butlerQC, inputRefs, outputRefs):
381 camera = butlerQC.get(inputRefs.camera)
382
383 if self.config.doOpticsTransmission:
384 opticsHandle = butlerQC.get(inputRefs.transmission_optics)
385 else:
386 opticsHandle = None
387
388 if self.config.doSensorTransmission:
389 sensorHandles = butlerQC.get(inputRefs.transmission_sensor)
390 sensorHandleDict = {sensorHandle.dataId['detector']: sensorHandle for
391 sensorHandle in sensorHandles}
392 else:
393 sensorHandleDict = {}
394
395 if self.config.doFilterDetectorTransmission:
396 filterHandles = butlerQC.get(inputRefs.transmission_filter_detector)
397 filterHandleDict = {
398 (filterHandle.dataId["physical_filter"], filterHandle.dataId["detector"]): filterHandle
399 for filterHandle in filterHandles
400 }
401 else:
402 filterHandles = butlerQC.get(inputRefs.transmission_filter)
403 filterHandleDict = {filterHandle.dataId['physical_filter']: filterHandle for
404 filterHandle in filterHandles}
405
406 filterToBand = {
407 filterHandle.dataId["physical_filter"]: filterHandle.dataId["band"]
408 for filterHandle in filterHandles
409 }
410
411 struct = self._fgcmMakeLut(
412 camera,
413 opticsHandle,
414 sensorHandleDict,
415 filterHandleDict,
416 filterToBand,
417 )
418
419 butlerQC.put(struct.fgcmLookUpTable, outputRefs.fgcmLookUpTable)
420 butlerQC.put(struct.fgcmStandardAtmosphere, outputRefs.fgcmStandardAtmosphere)
421
422 refDict = {passbandRef.dataId['physical_filter']: passbandRef for
423 passbandRef in outputRefs.fgcmStandardPassbands}
424 for physical_filter, passband in struct.fgcmStandardPassbands.items():
425 butlerQC.put(passband, refDict[physical_filter])
426
427 bandRefDict = {passbandRef.dataId["band"]: passbandRef for
428 passbandRef in outputRefs.standardPassbands}
429 for band, passband in struct.standardPassbands.items():
430 butlerQC.put(passband, bandRefDict[band])
431
432 def _fgcmMakeLut(self, camera, opticsHandle, sensorHandleDict,
433 filterHandleDict, filterToBand):
434 """
435 Make a FGCM Look-up Table
436
437 Parameters
438 ----------
439 camera : `lsst.afw.cameraGeom.Camera`
440 Camera from the butler.
441 opticsHandle : `lsst.daf.butler.DeferredDatasetHandle`
442 Reference to optics transmission curve.
443 sensorHandleDict : `dict` of [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
444 Dictionary of references to sensor transmission curves. Key will
445 be detector id.
446 filterHandleDict : `dict` of [`str`, `lsst.daf.butler.DeferredDatasetHandle`]
447 Dictionary of references to filter transmission curves. Key will
448 be physical filter label or tuple of physical filter label and
449 detector.
450 filterToBand : `dict` [`str`: `str`]
451 Mapping of physical filter to band name.
452
453 Returns
454 -------
455 retStruct : `lsst.pipe.base.Struct`
456 Output structure with keys:
457
458 fgcmLookUpTable : `BaseCatalog`
459 The FGCM look-up table.
460 fgcmStandardAtmosphere : `lsst.afw.image.TransmissionCurve`
461 Transmission curve for the FGCM standard atmosphere.
462 fgcmStandardPassbands : `dict` [`str`, `lsst.afw.image.TransmissionCurve`]
463 Dictionary of fgcm standard passbands, with the key as the
464 physical filter name.
465 standardPassbands : `dict` [`str`, `astropy.table.Table`]
466 Dictionary of standard passbands in astropy table format, with
467 the key as the band name.
468 """
469 # number of ccds from the length of the camera iterator
470 nCcd = countDetectors(camera, self.config.useScienceDetectors)
471
472 self.log.info("Found %d ccds for look-up table" % (nCcd))
473
474 # Load in optics, etc.
475 self._loadThroughputs(camera,
476 opticsHandle,
477 sensorHandleDict,
478 filterHandleDict)
479
480 lutConfig = self._createLutConfig(nCcd)
481
482 # make the lut object
483 self.log.info("Making the LUT maker object")
484 self.fgcmLutMaker = fgcm.FgcmLUTMaker(lutConfig)
485
486 # generate the throughput dictionary.
487
488 # these will be in Angstroms
489 # note that lambdaStep is currently in nm, because of historical
490 # reasons in the code. Convert to Angstroms here.
491 throughputLambda = np.arange(self.fgcmLutMaker.lambdaRange[0],
492 self.fgcmLutMaker.lambdaRange[1]+self.fgcmLutMaker.lambdaStep*10,
493 self.fgcmLutMaker.lambdaStep*10.)
494
495 self.log.info("Built throughput lambda, %.1f-%.1f, step %.2f" %
496 (throughputLambda[0], throughputLambda[-1],
497 throughputLambda[1] - throughputLambda[0]))
498
499 throughputDict = {}
500 for i, physicalFilter in enumerate(self.config.physicalFilters):
501 tDict = {}
502 tDict['LAMBDA'] = throughputLambda
503 for ccdIndex, detector in enumerate(camera):
504 if self.config.useScienceDetectors:
505 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE:
506 continue
507 tDict[ccdIndex] = self._getThroughputDetector(detector, physicalFilter, throughputLambda)
508 throughputDict[physicalFilter] = tDict
509
510 # set the throughputs
511 self.fgcmLutMaker.setThroughputs(throughputDict)
512
513 # make the LUT
514 self.log.info("Making LUT")
515 self.fgcmLutMaker.makeLUT()
516
517 # and save the LUT
518
519 # build the index values
520 comma = ','
521 physicalFilterString = comma.join(self.config.physicalFilters)
522 stdPhysicalFilterString = comma.join(self._getStdPhysicalFilterList())
523
524 atmosphereTableName = 'NoTableWasUsed'
525 if self.config.atmosphereTableName is not None:
526 atmosphereTableName = self.config.atmosphereTableName
527
528 lutSchema = self._makeLutSchema(physicalFilterString, stdPhysicalFilterString,
529 atmosphereTableName)
530
531 lutCat = self._makeLutCat(lutSchema, physicalFilterString,
532 stdPhysicalFilterString, atmosphereTableName)
533
534 atmStd = TransmissionCurve.makeSpatiallyConstant(
535 throughput=self.fgcmLutMaker.atmStdTrans.astype(np.float64),
536 wavelengths=self.fgcmLutMaker.atmLambda.astype(np.float64),
537 throughputAtMin=self.fgcmLutMaker.atmStdTrans[0],
538 throughputAtMax=self.fgcmLutMaker.atmStdTrans[1],
539 )
540
541 fgcmStandardPassbands = {}
542 for i, physical_filter in enumerate(self.fgcmLutMaker.filterNames):
543 passband = self.fgcmLutMaker.throughputs[i]['THROUGHPUT_AVG']*self.fgcmLutMaker.atmStdTrans
544 fgcmStandardPassbands[physical_filter] = TransmissionCurve.makeSpatiallyConstant(
545 throughput=passband.astype(np.float64),
546 wavelengths=self.fgcmLutMaker.atmLambda.astype(np.float64),
547 throughputAtMin=passband[0],
548 throughputAtMax=passband[-1],
549 )
550
551 standardPassbands = {}
552 for i, physical_filter in enumerate(self.fgcmLutMaker.filterNames):
553 if physical_filter != self.fgcmLutMaker.stdFilterNames[i]:
554 # This filter does not map onto one of the "standard"
555 # passbands. E.g., for HSC we have HSC-R and HSC-R2 which
556 # both map onto HSC-R2 which sets the "standard".
557 continue
558 band = filterToBand[physical_filter]
559 passband = self.fgcmLutMaker.throughputs[i]['THROUGHPUT_AVG']*self.fgcmLutMaker.atmStdTrans
560 passbandTable = Table(
561 {
562 "wavelength": (self.fgcmLutMaker.atmLambda.astype(np.float64)/10.)*units.nm,
563 "throughput": (passband*100.)*units.percent,
564 },
565 )
566 passbandTable["wavelength"].description = "Wavelength bin centers"
567 standardPassbands[band] = passbandTable
568
569 retStruct = pipeBase.Struct(
570 fgcmLookUpTable=lutCat,
571 fgcmStandardAtmosphere=atmStd,
572 fgcmStandardPassbands=fgcmStandardPassbands,
573 standardPassbands=standardPassbands,
574 )
575
576 return retStruct
577
579 """Get the standard physical filter lists from config.physicalFilters
580 and config.stdPhysicalFilterOverrideMap
581
582 Returns
583 -------
584 stdPhysicalFilters : `list`
585 """
586 override = self.config.stdPhysicalFilterOverrideMap
587 return [override.get(physicalFilter, physicalFilter) for
588 physicalFilter in self.config.physicalFilters]
589
590 def _createLutConfig(self, nCcd):
591 """
592 Create the fgcmLut config dictionary
593
594 Parameters
595 ----------
596 nCcd: `int`
597 Number of CCDs in the camera
598 """
599
600 # create the common stub of the lutConfig
601 lutConfig = {}
602 lutConfig['logger'] = self.log
603 lutConfig['filterNames'] = self.config.physicalFilters
604 lutConfig['stdFilterNames'] = self._getStdPhysicalFilterList()
605 lutConfig['nCCD'] = nCcd
606
607 # atmosphereTable already validated if available
608 if self.config.atmosphereTableName is not None:
609 lutConfig['atmosphereTableName'] = self.config.atmosphereTableName
610 else:
611 # use the regular paramters (also validated if needed)
612 lutConfig['elevation'] = self.config.parameters.elevation
613 lutConfig['pmbRange'] = self.config.parameters.pmbRange
614 lutConfig['pmbSteps'] = self.config.parameters.pmbSteps
615 lutConfig['pwvRange'] = self.config.parameters.pwvRange
616 lutConfig['pwvSteps'] = self.config.parameters.pwvSteps
617 lutConfig['o3Range'] = self.config.parameters.o3Range
618 lutConfig['o3Steps'] = self.config.parameters.o3Steps
619 lutConfig['tauRange'] = self.config.parameters.tauRange
620 lutConfig['tauSteps'] = self.config.parameters.tauSteps
621 lutConfig['alphaRange'] = self.config.parameters.alphaRange
622 lutConfig['alphaSteps'] = self.config.parameters.alphaSteps
623 lutConfig['zenithRange'] = self.config.parameters.zenithRange
624 lutConfig['zenithSteps'] = self.config.parameters.zenithSteps
625 lutConfig['pmbStd'] = self.config.parameters.pmbStd
626 lutConfig['pwvStd'] = self.config.parameters.pwvStd
627 lutConfig['o3Std'] = self.config.parameters.o3Std
628 lutConfig['tauStd'] = self.config.parameters.tauStd
629 lutConfig['alphaStd'] = self.config.parameters.alphaStd
630 lutConfig['airmassStd'] = self.config.parameters.airmassStd
631 lutConfig['lambdaRange'] = self.config.parameters.lambdaRange
632 lutConfig['lambdaStep'] = self.config.parameters.lambdaStep
633 lutConfig['lambdaNorm'] = self.config.parameters.lambdaNorm
634
635 # Add any per-filter correction term updates if necessary.
636 # Note that sensorCTerms is the name of the config field in fgcm.
637 if self.config.sensorCorrectionTermDict:
638 lutConfig['sensorCTerms'] = {}
639 for key, value in self.config.sensorCorrectionTermDict.items():
640 lutConfig['sensorCTerms'][key] = (
641 value.refLambda,
642 dict(value.correctionTermDict),
643 )
644
645 return lutConfig
646
647 def _loadThroughputs(self, camera, opticsHandle, sensorHandleDict, filterHandleDict):
648 """Internal method to load throughput data for filters
649
650 Parameters
651 ----------
652 camera: `lsst.afw.cameraGeom.Camera`
653 Camera from the butler
654 opticsHandle : `lsst.daf.butler.DeferredDatasetHandle`
655 Reference to optics transmission curve.
656 sensorHandleDict : `dict` of [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
657 Dictionary of references to sensor transmission curves. Key will
658 be detector id.
659 filterHandleDict : `dict` of [`str`, `lsst.daf.butler.DeferredDatasetHandle`]
660 Dictionary of references to filter transmission curves. Key will
661 be physical filter label.
662
663 Raises
664 ------
665 ValueError : Raised if configured filter name does not match any of the
666 available filter transmission curves.
667 """
668 if self.config.doOpticsTransmission:
669 self._opticsTransmission = opticsHandle.get()
670 else:
671 self._opticsTransmission = TransmissionCurve.makeSpatiallyConstant(
672 throughput=np.ones(100),
673 wavelengths=np.linspace(
674 self.config.parameters.lambdaRange[0],
675 self.config.parameters.lambdaRange[1],
676 100,
677 ),
678 throughputAtMin=1.0,
679 throughputAtMax=1.0,
680 )
681
683 for detector in camera:
684 if self.config.useScienceDetectors:
685 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE:
686 continue
687 if self.config.doSensorTransmission:
688 self._sensorsTransmission[detector.getId()] = sensorHandleDict[detector.getId()].get()
689 else:
690 self._sensorsTransmission[detector.getId()] = TransmissionCurve.makeSpatiallyConstant(
691 throughput=np.ones(100),
692 wavelengths=np.linspace(
693 self.config.parameters.lambdaRange[0],
694 self.config.parameters.lambdaRange[1],
695 100,
696 ),
697 throughputAtMin=1.0,
698 throughputAtMax=1.0,
699 )
700
702 if self.config.doFilterDetectorTransmission:
703 for physicalFilter in self.config.physicalFilters:
704 for detector in camera:
705 if self.config.useScienceDetectors:
706 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE:
707 continue
708 key = (physicalFilter, detector.getId())
709 self._filtersTransmission[key] = filterHandleDict[key].get()
710 else:
711 for physicalFilter in self.config.physicalFilters:
712 self._filtersTransmission[physicalFilter] = filterHandleDict[physicalFilter].get()
713
714 def _getThroughputDetector(self, detector, physicalFilter, throughputLambda):
715 """Internal method to get throughput for a detector.
716
717 Returns the throughput at the center of the detector for a given filter.
718
719 Parameters
720 ----------
721 detector: `lsst.afw.cameraGeom._detector.Detector`
722 Detector on camera
723 physicalFilter: `str`
724 Physical filter label
725 throughputLambda: `np.array(dtype=np.float64)`
726 Wavelength steps (Angstrom)
727
728 Returns
729 -------
730 throughput: `np.array(dtype=np.float64)`
731 Throughput (max 1.0) at throughputLambda
732 """
733
734 c = detector.getCenter(afwCameraGeom.FOCAL_PLANE)
735 c.scale(1.0/detector.getPixelSize()[0]) # Assumes x and y pixel sizes in arcsec are the same
736
737 throughput = self._opticsTransmission.sampleAt(position=c,
738 wavelengths=throughputLambda)
739
740 throughput *= self._sensorsTransmission[detector.getId()].sampleAt(position=c,
741 wavelengths=throughputLambda)
742
743 if self.config.doFilterDetectorTransmission:
744 throughput *= self._filtersTransmission[(physicalFilter, detector.getId())].sampleAt(
745 position=c,
746 wavelengths=throughputLambda,
747 )
748 else:
749 throughput *= self._filtersTransmission[physicalFilter].sampleAt(
750 position=c,
751 wavelengths=throughputLambda,
752 )
753
754 # Clip the throughput from 0 to 1
755 throughput = np.clip(throughput, 0.0, 1.0)
756
757 return throughput
758
759 def _makeLutSchema(self, physicalFilterString, stdPhysicalFilterString,
760 atmosphereTableName):
761 """
762 Make the LUT schema
763
764 Parameters
765 ----------
766 physicalFilterString: `str`
767 Combined string of all the physicalFilters
768 stdPhysicalFilterString: `str`
769 Combined string of all the standard physicalFilters
770 atmosphereTableName: `str`
771 Name of the atmosphere table used to generate LUT
772
773 Returns
774 -------
775 lutSchema: `afwTable.schema`
776 """
777
778 lutSchema = afwTable.Schema()
779
780 lutSchema.addField('tablename', type=str, doc='Atmosphere table name',
781 size=len(atmosphereTableName))
782 lutSchema.addField('elevation', type=float, doc="Telescope elevation used for LUT")
783 lutSchema.addField('physicalFilters', type=str, doc='physicalFilters in LUT',
784 size=len(physicalFilterString))
785 lutSchema.addField('stdPhysicalFilters', type=str, doc='Standard physicalFilters in LUT',
786 size=len(stdPhysicalFilterString))
787 lutSchema.addField('pmb', type='ArrayD', doc='Barometric Pressure',
788 size=self.fgcmLutMaker.pmb.size)
789 lutSchema.addField('pmbFactor', type='ArrayD', doc='PMB scaling factor',
790 size=self.fgcmLutMaker.pmb.size)
791 lutSchema.addField('pmbElevation', type=np.float64, doc='PMB Scaling at elevation')
792 lutSchema.addField('pwv', type='ArrayD', doc='Preciptable Water Vapor',
793 size=self.fgcmLutMaker.pwv.size)
794 lutSchema.addField('o3', type='ArrayD', doc='Ozone',
795 size=self.fgcmLutMaker.o3.size)
796 lutSchema.addField('tau', type='ArrayD', doc='Aerosol optical depth',
797 size=self.fgcmLutMaker.tau.size)
798 lutSchema.addField('lambdaNorm', type=np.float64, doc='AOD wavelength')
799 lutSchema.addField('alpha', type='ArrayD', doc='Aerosol alpha',
800 size=self.fgcmLutMaker.alpha.size)
801 lutSchema.addField('zenith', type='ArrayD', doc='Zenith angle',
802 size=self.fgcmLutMaker.zenith.size)
803 lutSchema.addField('nCcd', type=np.int32, doc='Number of CCDs')
804
805 # and the standard values
806 lutSchema.addField('pmbStd', type=np.float64, doc='PMB Standard')
807 lutSchema.addField('pwvStd', type=np.float64, doc='PWV Standard')
808 lutSchema.addField('o3Std', type=np.float64, doc='O3 Standard')
809 lutSchema.addField('tauStd', type=np.float64, doc='Tau Standard')
810 lutSchema.addField('alphaStd', type=np.float64, doc='Alpha Standard')
811 lutSchema.addField('zenithStd', type=np.float64, doc='Zenith angle Standard')
812 lutSchema.addField('lambdaRange', type='ArrayD', doc='Wavelength range',
813 size=2)
814 lutSchema.addField('lambdaStep', type=np.float64, doc='Wavelength step')
815 lutSchema.addField('lambdaStd', type='ArrayD', doc='Standard Wavelength',
816 size=len(self.fgcmLutMaker.filterNames))
817 lutSchema.addField('lambdaStdFilter', type='ArrayD', doc='Standard Wavelength (raw)',
818 size=len(self.fgcmLutMaker.filterNames))
819 lutSchema.addField('i0Std', type='ArrayD', doc='I0 Standard',
820 size=len(self.fgcmLutMaker.filterNames))
821 lutSchema.addField('i1Std', type='ArrayD', doc='I1 Standard',
822 size=len(self.fgcmLutMaker.filterNames))
823 lutSchema.addField('i10Std', type='ArrayD', doc='I10 Standard',
824 size=len(self.fgcmLutMaker.filterNames))
825 lutSchema.addField('i2Std', type='ArrayD', doc='I2 Standard',
826 size=len(self.fgcmLutMaker.filterNames))
827 lutSchema.addField('lambdaB', type='ArrayD', doc='Wavelength for passband (no atm)',
828 size=len(self.fgcmLutMaker.filterNames))
829 lutSchema.addField('atmLambda', type='ArrayD', doc='Atmosphere wavelengths (Angstrom)',
830 size=self.fgcmLutMaker.atmLambda.size)
831 lutSchema.addField('atmStdTrans', type='ArrayD', doc='Standard Atmosphere Throughput',
832 size=self.fgcmLutMaker.atmStdTrans.size)
833
834 # and the look-up-tables
835 lutSchema.addField('luttype', type=str, size=20, doc='Look-up table type')
836 lutSchema.addField('lut', type='ArrayF', doc='Look-up table for luttype',
837 size=self.fgcmLutMaker.lut['I0'].size)
838
839 return lutSchema
840
841 def _makeLutCat(self, lutSchema, physicalFilterString, stdPhysicalFilterString,
842 atmosphereTableName):
843 """
844 Make the LUT schema
845
846 Parameters
847 ----------
848 lutSchema: `afwTable.schema`
849 Lut catalog schema
850 physicalFilterString: `str`
851 Combined string of all the physicalFilters
852 stdPhysicalFilterString: `str`
853 Combined string of all the standard physicalFilters
854 atmosphereTableName: `str`
855 Name of the atmosphere table used to generate LUT
856
857 Returns
858 -------
859 lutCat: `afwTable.BaseCatalog`
860 Look-up table catalog for persistence.
861 """
862
863 # The somewhat strange format is to make sure that
864 # the rows of the afwTable do not get too large
865 # (see DM-11419)
866
867 lutCat = afwTable.BaseCatalog(lutSchema)
868 lutCat.table.preallocate(14)
869
870 # first fill the first index
871 rec = lutCat.addNew()
872
873 rec['tablename'] = atmosphereTableName
874 rec['elevation'] = self.fgcmLutMaker.atmosphereTable.elevation
875 rec['physicalFilters'] = physicalFilterString
876 rec['stdPhysicalFilters'] = stdPhysicalFilterString
877 rec['pmb'][:] = self.fgcmLutMaker.pmb
878 rec['pmbFactor'][:] = self.fgcmLutMaker.pmbFactor
879 rec['pmbElevation'] = self.fgcmLutMaker.pmbElevation
880 rec['pwv'][:] = self.fgcmLutMaker.pwv
881 rec['o3'][:] = self.fgcmLutMaker.o3
882 rec['tau'][:] = self.fgcmLutMaker.tau
883 rec['lambdaNorm'] = self.fgcmLutMaker.lambdaNorm
884 rec['alpha'][:] = self.fgcmLutMaker.alpha
885 rec['zenith'][:] = self.fgcmLutMaker.zenith
886 rec['nCcd'] = self.fgcmLutMaker.nCCD
887
888 rec['pmbStd'] = self.fgcmLutMaker.pmbStd
889 rec['pwvStd'] = self.fgcmLutMaker.pwvStd
890 rec['o3Std'] = self.fgcmLutMaker.o3Std
891 rec['tauStd'] = self.fgcmLutMaker.tauStd
892 rec['alphaStd'] = self.fgcmLutMaker.alphaStd
893 rec['zenithStd'] = self.fgcmLutMaker.zenithStd
894 rec['lambdaRange'][:] = self.fgcmLutMaker.lambdaRange
895 rec['lambdaStep'] = self.fgcmLutMaker.lambdaStep
896 rec['lambdaStd'][:] = self.fgcmLutMaker.lambdaStd
897 rec['lambdaStdFilter'][:] = self.fgcmLutMaker.lambdaStdFilter
898 rec['i0Std'][:] = self.fgcmLutMaker.I0Std
899 rec['i1Std'][:] = self.fgcmLutMaker.I1Std
900 rec['i10Std'][:] = self.fgcmLutMaker.I10Std
901 rec['i2Std'][:] = self.fgcmLutMaker.I2Std
902 rec['lambdaB'][:] = self.fgcmLutMaker.lambdaB
903 rec['atmLambda'][:] = self.fgcmLutMaker.atmLambda
904 rec['atmStdTrans'][:] = self.fgcmLutMaker.atmStdTrans
905
906 rec['luttype'] = 'I0'
907 rec['lut'][:] = self.fgcmLutMaker.lut['I0'].flatten()
908
909 # and add the rest
910 rec = lutCat.addNew()
911 rec['luttype'] = 'I1'
912 rec['lut'][:] = self.fgcmLutMaker.lut['I1'].flatten()
913
914 derivTypes = ['D_PMB', 'D_LNPWV', 'D_O3', 'D_LNTAU', 'D_ALPHA', 'D_SECZENITH',
915 'D_PMB_I1', 'D_LNPWV_I1', 'D_O3_I1', 'D_LNTAU_I1', 'D_ALPHA_I1',
916 'D_SECZENITH_I1']
917 for derivType in derivTypes:
918 rec = lutCat.addNew()
919 rec['luttype'] = derivType
920 rec['lut'][:] = self.fgcmLutMaker.lutDeriv[derivType].flatten()
921
922 return lutCat
Defines the fields and offsets for a table.
Definition Schema.h:51
_makeLutSchema(self, physicalFilterString, stdPhysicalFilterString, atmosphereTableName)
_getThroughputDetector(self, detector, physicalFilter, throughputLambda)
_loadThroughputs(self, camera, opticsHandle, sensorHandleDict, filterHandleDict)
_fgcmMakeLut(self, camera, opticsHandle, sensorHandleDict, filterHandleDict, filterToBand)
runQuantum(self, butlerQC, inputRefs, outputRefs)
__init__(self, initInputs=None, **kwargs)
_makeLutCat(self, lutSchema, physicalFilterString, stdPhysicalFilterString, atmosphereTableName)