LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
fgcmFitCycle.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"""Perform a single fit cycle of FGCM.
24
25This task runs a single "fit cycle" of fgcm. Prior to running this task
26one must run both fgcmMakeLut (to construct the atmosphere and instrumental
27look-up-table) and fgcmBuildStars (to extract visits and star observations
28for the global fit).
29
30The fgcmFitCycle is meant to be run multiple times, and is tracked by the
31'cycleNumber'. After each run of the fit cycle, diagnostic plots should
32be inspected to set parameters for outlier rejection on the following
33cycle. Please see the fgcmcal Cookbook for details.
34"""
35
36import copy
37
38import numpy as np
39
40import lsst.pex.config as pexConfig
41import lsst.pipe.base as pipeBase
42from lsst.pipe.base import connectionTypes
43import lsst.afw.table as afwTable
44
45from .utilities import makeConfigDict, translateFgcmLut, translateVisitCatalog
46from .utilities import extractReferenceMags
47from .utilities import makeZptSchema, makeZptCat
48from .utilities import makeAtmSchema, makeAtmCat, makeStdSchema, makeStdCat
49from .sedterms import SedboundarytermDict, SedtermDict
50from .focalPlaneProjector import FocalPlaneProjector
51
52import fgcm
53
54__all__ = ['FgcmFitCycleConfig', 'FgcmFitCycleTask']
55
56MULTIPLE_CYCLES_MAX = 10
57
58
59class FgcmFitCycleConnections(pipeBase.PipelineTaskConnections,
60 dimensions=("instrument",),
61 defaultTemplates={"previousCycleNumber": "-1",
62 "cycleNumber": "0"}):
63 camera = connectionTypes.PrerequisiteInput(
64 doc="Camera instrument",
65 name="camera",
66 storageClass="Camera",
67 dimensions=("instrument",),
68 isCalibration=True,
69 )
70
71 fgcmLookUpTable = connectionTypes.PrerequisiteInput(
72 doc=("Atmosphere + instrument look-up-table for FGCM throughput and "
73 "chromatic corrections."),
74 name="fgcmLookUpTable",
75 storageClass="Catalog",
76 dimensions=("instrument",),
77 deferLoad=True,
78 )
79
80 fgcmVisitCatalog = connectionTypes.Input(
81 doc="Catalog of visit information for fgcm",
82 name="fgcmVisitCatalog",
83 storageClass="Catalog",
84 dimensions=("instrument",),
85 deferLoad=True,
86 )
87
88 fgcmStarObservationsParquet = connectionTypes.Input(
89 doc=("Catalog of star observations for fgcm, in parquet format. "
90 "Used if useParquetCatalogFormat is True."),
91 name="fgcm_star_observations",
92 storageClass="ArrowAstropy",
93 dimensions=("instrument",),
94 deferLoad=True,
95 )
96
97 fgcmStarIdsParquet = connectionTypes.Input(
98 doc=("Catalog of fgcm calibration star IDs, in parquet format. "
99 "Used if useParquetCatalogFormat is True."),
100 name="fgcm_star_ids",
101 storageClass="ArrowAstropy",
102 dimensions=("instrument",),
103 deferLoad=True,
104 )
105
106 fgcmReferenceStarsParquet = connectionTypes.Input(
107 doc=("Catalog of fgcm-matched reference stars, in parquet format. "
108 "Used if useParquetCatalogFormat is True."),
109 name="fgcm_reference_stars",
110 storageClass="ArrowAstropy",
111 dimensions=("instrument",),
112 deferLoad=True,
113 )
114
115 fgcmStarObservations = connectionTypes.Input(
116 doc=("Catalog of star observations for fgcm; old format. "
117 "Used if useParquetCatalogFormat is False."),
118 name="fgcmStarObservations",
119 storageClass="Catalog",
120 dimensions=("instrument",),
121 deferLoad=True,
122 )
123
124 fgcmStarIds = connectionTypes.Input(
125 doc=("Catalog of fgcm calibration star IDs. "
126 "Used if useParquetCatalogFormat is False."),
127 name="fgcmStarIds",
128 storageClass="Catalog",
129 dimensions=("instrument",),
130 deferLoad=True,
131 )
132
133 fgcmStarIndices = connectionTypes.Input(
134 doc=("Catalog of fgcm calibration star indices; old format."
135 "Used if useParquetCatalogFormat is False."),
136 name="fgcmStarIndices",
137 storageClass="Catalog",
138 dimensions=("instrument",),
139 deferLoad=True,
140 )
141
142 fgcmReferenceStars = connectionTypes.Input(
143 doc=("Catalog of fgcm-matched reference stars; old format."
144 "Used if useParquetCatalogFormat is False."),
145 name="fgcmReferenceStars",
146 storageClass="Catalog",
147 dimensions=("instrument",),
148 deferLoad=True,
149 )
150
151 fgcmFlaggedStarsInput = connectionTypes.PrerequisiteInput(
152 doc="Catalog of flagged stars for fgcm calibration from previous fit cycle",
153 name="fgcmFlaggedStars{previousCycleNumber}",
154 storageClass="Catalog",
155 dimensions=("instrument",),
156 deferLoad=True,
157 )
158
159 fgcmFitParametersInput = connectionTypes.PrerequisiteInput(
160 doc="Catalog of fgcm fit parameters from previous fit cycle",
161 name="fgcmFitParameters{previousCycleNumber}",
162 storageClass="Catalog",
163 dimensions=("instrument",),
164 deferLoad=True,
165 )
166
167 fgcmFitParameters = connectionTypes.Output(
168 doc="Catalog of fgcm fit parameters from current fit cycle",
169 name="fgcmFitParameters{cycleNumber}",
170 storageClass="Catalog",
171 dimensions=("instrument",),
172 )
173
174 fgcmFlaggedStars = connectionTypes.Output(
175 doc="Catalog of flagged stars for fgcm calibration from current fit cycle",
176 name="fgcmFlaggedStars{cycleNumber}",
177 storageClass="Catalog",
178 dimensions=("instrument",),
179 )
180
181 fgcmZeropoints = connectionTypes.Output(
182 doc="Catalog of fgcm zeropoint data from current fit cycle",
183 name="fgcmZeropoints{cycleNumber}",
184 storageClass="Catalog",
185 dimensions=("instrument",),
186 )
187
188 fgcmAtmosphereParameters = connectionTypes.Output(
189 doc="Catalog of atmospheric fit parameters from current fit cycle",
190 name="fgcmAtmosphereParameters{cycleNumber}",
191 storageClass="Catalog",
192 dimensions=("instrument",),
193 )
194
195 fgcmStandardStars = connectionTypes.Output(
196 doc="Catalog of standard star magnitudes from current fit cycle",
197 name="fgcmStandardStars{cycleNumber}",
198 storageClass="SimpleCatalog",
199 dimensions=("instrument",),
200 )
201
202 # Add connections for running multiple cycles
203 # This uses vars() to alter the class __dict__ to programmatically
204 # write many similar outputs.
205 for cycle in range(MULTIPLE_CYCLES_MAX):
206 vars()[f"fgcmFitParameters{cycle}"] = connectionTypes.Output(
207 doc=f"Catalog of fgcm fit parameters from fit cycle {cycle}",
208 name=f"fgcmFitParameters{cycle}",
209 storageClass="Catalog",
210 dimensions=("instrument",),
211 )
212 vars()[f"fgcmFlaggedStars{cycle}"] = connectionTypes.Output(
213 doc=f"Catalog of flagged stars for fgcm calibration from fit cycle {cycle}",
214 name=f"fgcmFlaggedStars{cycle}",
215 storageClass="Catalog",
216 dimensions=("instrument",),
217 )
218 vars()[f"fgcmZeropoints{cycle}"] = connectionTypes.Output(
219 doc=f"Catalog of fgcm zeropoint data from fit cycle {cycle}",
220 name=f"fgcmZeropoints{cycle}",
221 storageClass="Catalog",
222 dimensions=("instrument",),
223 )
224 vars()[f"fgcmAtmosphereParameters{cycle}"] = connectionTypes.Output(
225 doc=f"Catalog of atmospheric fit parameters from fit cycle {cycle}",
226 name=f"fgcmAtmosphereParameters{cycle}",
227 storageClass="Catalog",
228 dimensions=("instrument",),
229 )
230 vars()[f"fgcmStandardStars{cycle}"] = connectionTypes.Output(
231 doc=f"Catalog of standard star magnitudes from fit cycle {cycle}",
232 name=f"fgcmStandardStars{cycle}",
233 storageClass="SimpleCatalog",
234 dimensions=("instrument",),
235 )
236
237 def __init__(self, *, config=None):
238 super().__init__(config=config)
239
240 if not config.doReferenceCalibration:
241 self.inputs.remove("fgcmReferenceStars")
242 self.inputs.remove("fgcmReferenceStarsParquet")
243
244 if config.useParquetCatalogFormat:
245 self.inputs.remove("fgcmStarObservations")
246 self.inputs.remove("fgcmStarIds")
247 self.inputs.remove("fgcmStarIndices")
248 if config.doReferenceCalibration:
249 self.inputs.remove("fgcmReferenceStars")
250 else:
251 self.inputs.remove("fgcmStarObservationsParquet")
252 self.inputs.remove("fgcmStarIdsParquet")
253 if config.doReferenceCalibration:
254 self.inputs.remove("fgcmReferenceStarsParquet")
255
256 if str(int(config.connections.cycleNumber)) != config.connections.cycleNumber:
257 raise ValueError("cycleNumber must be of integer format")
258 if str(int(config.connections.previousCycleNumber)) != config.connections.previousCycleNumber:
259 raise ValueError("previousCycleNumber must be of integer format")
260 if int(config.connections.previousCycleNumber) != (int(config.connections.cycleNumber) - 1):
261 raise ValueError("previousCycleNumber must be 1 less than cycleNumber")
262
263 if int(config.connections.cycleNumber) == 0:
264 self.prerequisiteInputs.remove("fgcmFlaggedStarsInput")
265 self.prerequisiteInputs.remove("fgcmFitParametersInput")
266
267 if not self.config.doMultipleCycles:
268 # Single-cycle run
269 if not self.config.isFinalCycle and not self.config.outputStandardsBeforeFinalCycle:
270 self.outputs.remove("fgcmStandardStars")
271
272 if not self.config.isFinalCycle and not self.config.outputZeropointsBeforeFinalCycle:
273 self.outputs.remove("fgcmZeropoints")
274 self.outputs.remove("fgcmAtmosphereParameters")
275
276 # Remove all multiple cycle outputs
277 for cycle in range(0, MULTIPLE_CYCLES_MAX):
278 self.outputs.remove(f"fgcmFitParameters{cycle}")
279 self.outputs.remove(f"fgcmFlaggedStars{cycle}")
280 self.outputs.remove(f"fgcmZeropoints{cycle}")
281 self.outputs.remove(f"fgcmAtmosphereParameters{cycle}")
282 self.outputs.remove(f"fgcmStandardStars{cycle}")
283
284 else:
285 # Multiple-cycle run
286 # Remove single-cycle outputs
287 self.outputs.remove("fgcmFitParameters")
288 self.outputs.remove("fgcmFlaggedStars")
289 self.outputs.remove("fgcmZeropoints")
290 self.outputs.remove("fgcmAtmosphereParameters")
291 self.outputs.remove("fgcmStandardStars")
292
293 # Remove outputs from cycles that are not used
294 for cycle in range(self.config.multipleCyclesFinalCycleNumber + 1,
295 MULTIPLE_CYCLES_MAX):
296 self.outputs.remove(f"fgcmFitParameters{cycle}")
297 self.outputs.remove(f"fgcmFlaggedStars{cycle}")
298 self.outputs.remove(f"fgcmZeropoints{cycle}")
299 self.outputs.remove(f"fgcmAtmosphereParameters{cycle}")
300 self.outputs.remove(f"fgcmStandardStars{cycle}")
301
302 # Remove non-final-cycle outputs if necessary
303 for cycle in range(self.config.multipleCyclesFinalCycleNumber):
304 if not self.config.outputZeropointsBeforeFinalCycle:
305 self.outputs.remove(f"fgcmZeropoints{cycle}")
306 self.outputs.remove(f"fgcmAtmosphereParameters{cycle}")
307 if not self.config.outputStandardsBeforeFinalCycle:
308 self.outputs.remove(f"fgcmStandardStars{cycle}")
309
310
311class FgcmFitCycleConfig(pipeBase.PipelineTaskConfig,
312 pipelineConnections=FgcmFitCycleConnections):
313 """Config for FgcmFitCycle"""
314
315 doMultipleCycles = pexConfig.Field(
316 doc="Run multiple fit cycles in one task",
317 dtype=bool,
318 default=False,
319 )
320 useParquetCatalogFormat = pexConfig.Field(
321 doc="Use parquet catalog format?",
322 dtype=bool,
323 default=True,
324 )
325 multipleCyclesFinalCycleNumber = pexConfig.RangeField(
326 doc=("Final cycle number in multiple cycle mode. The initial cycle "
327 "is 0, with limited parameters fit. The next cycle is 1 with "
328 "full parameter fit. The final cycle is a clean-up with no "
329 "parameters fit. There will be a total of "
330 "(multipleCycleFinalCycleNumber + 1) cycles run, and the final "
331 "cycle number cannot be less than 2."),
332 dtype=int,
333 default=5,
334 min=2,
335 max=MULTIPLE_CYCLES_MAX,
336 inclusiveMax=True,
337 )
338 bands = pexConfig.ListField(
339 doc="Bands to run calibration",
340 dtype=str,
341 default=[],
342 )
343 fitBands = pexConfig.ListField(
344 doc=("Bands to use in atmospheric fit. The bands not listed here will have "
345 "the atmosphere constrained from the 'fitBands' on the same night. "
346 "Must be a subset of `config.bands`"),
347 dtype=str,
348 default=[],
349 )
350 requiredBands = pexConfig.ListField(
351 doc=("Bands that are required for a star to be considered a calibration star. "
352 "Must be a subset of `config.bands`"),
353 dtype=str,
354 default=[],
355 )
356 physicalFilterMap = pexConfig.DictField(
357 doc="Mapping from 'physicalFilter' to band.",
358 keytype=str,
359 itemtype=str,
360 default={},
361 )
362 doReferenceCalibration = pexConfig.Field(
363 doc="Use reference catalog as additional constraint on calibration",
364 dtype=bool,
365 default=True,
366 )
367 refStarSnMin = pexConfig.Field(
368 doc="Reference star signal-to-noise minimum to use in calibration. Set to <=0 for no cut.",
369 dtype=float,
370 default=50.0,
371 )
372 refStarOutlierNSig = pexConfig.Field(
373 doc=("Number of sigma compared to average mag for reference star to be considered an outlier. "
374 "Computed per-band, and if it is an outlier in any band it is rejected from fits."),
375 dtype=float,
376 default=4.0,
377 )
378 applyRefStarColorCuts = pexConfig.Field(
379 doc=("Apply color cuts defined in ``starColorCuts`` to reference stars? "
380 "These cuts are in addition to any cuts defined in ``refStarColorCuts``"),
381 dtype=bool,
382 default=True,
383 )
384 refStarMaxFracUse = pexConfig.Field(
385 doc=("Maximum fraction of reference stars to use in the fit. Remainder will "
386 "be used only for validation."),
387 dtype=float,
388 default=0.5,
389 )
390 useExposureReferenceOffset = pexConfig.Field(
391 doc=("Use per-exposure (visit) offsets between calibrated stars and reference stars "
392 "for final zeropoints? This may help uniformity for disjoint surveys."),
393 dtype=bool,
394 default=False,
395 )
396 nCore = pexConfig.Field(
397 doc="Number of cores to use",
398 dtype=int,
399 default=4,
400 deprecated="Number of cores is deprecated as a config, and will be removed after v27. "
401 "Please use ``pipetask run --cores-per-quantum`` instead.",
402 )
403 nStarPerRun = pexConfig.Field(
404 doc="Number of stars to run in each chunk",
405 dtype=int,
406 default=200000,
407 )
408 nExpPerRun = pexConfig.Field(
409 doc="Number of exposures to run in each chunk",
410 dtype=int,
411 default=1000,
412 )
413 reserveFraction = pexConfig.Field(
414 doc="Fraction of stars to reserve for testing",
415 dtype=float,
416 default=0.1,
417 )
418 freezeStdAtmosphere = pexConfig.Field(
419 doc="Freeze atmosphere parameters to standard (for testing)",
420 dtype=bool,
421 default=False,
422 )
423 precomputeSuperStarInitialCycle = pexConfig.Field(
424 doc="Precompute superstar flat for initial cycle",
425 dtype=bool,
426 default=False,
427 )
428 superStarSubCcdDict = pexConfig.DictField(
429 doc=("Per-band specification on whether to compute superstar flat on sub-ccd scale. "
430 "Must have one entry per band."),
431 keytype=str,
432 itemtype=bool,
433 default={},
434 )
435 superStarSubCcdChebyshevOrder = pexConfig.Field(
436 doc=("Order of the 2D chebyshev polynomials for sub-ccd superstar fit. "
437 "Global default is first-order polynomials, and should be overridden "
438 "on a camera-by-camera basis depending on the ISR."),
439 dtype=int,
440 default=1,
441 )
442 superStarSubCcdTriangular = pexConfig.Field(
443 doc=("Should the sub-ccd superstar chebyshev matrix be triangular to "
444 "suppress high-order cross terms?"),
445 dtype=bool,
446 default=False,
447 )
448 superStarSigmaClip = pexConfig.Field(
449 doc="Number of sigma to clip outliers when selecting for superstar flats",
450 dtype=float,
451 default=5.0,
452 )
453 superStarPlotCcdResiduals = pexConfig.Field(
454 doc="If plotting is enabled, should per-detector residuals be plotted? "
455 "This may produce a lot of output, and should be used only for "
456 "debugging purposes.",
457 dtype=bool,
458 default=False,
459 )
460 focalPlaneSigmaClip = pexConfig.Field(
461 doc="Number of sigma to clip outliers per focal-plane.",
462 dtype=float,
463 default=4.0,
464 )
465 ccdGraySubCcdDict = pexConfig.DictField(
466 doc=("Per-band specification on whether to compute achromatic per-ccd residual "
467 "('ccd gray') on a sub-ccd scale."),
468 keytype=str,
469 itemtype=bool,
470 default={},
471 )
472 ccdGraySubCcdChebyshevOrder = pexConfig.Field(
473 doc="Order of the 2D chebyshev polynomials for sub-ccd gray fit.",
474 dtype=int,
475 default=1,
476 )
477 ccdGraySubCcdTriangular = pexConfig.Field(
478 doc=("Should the sub-ccd gray chebyshev matrix be triangular to "
479 "suppress high-order cross terms?"),
480 dtype=bool,
481 default=True,
482 )
483 ccdGrayFocalPlaneDict = pexConfig.DictField(
484 doc=("Per-band specification on whether to compute focal-plane residual "
485 "('ccd gray') corrections."),
486 keytype=str,
487 itemtype=bool,
488 default={},
489 )
490 ccdGrayFocalPlaneFitMinCcd = pexConfig.Field(
491 doc=("Minimum number of 'good' CCDs required to perform focal-plane "
492 "gray corrections. If there are fewer good CCDs then the gray "
493 "correction is computed per-ccd."),
494 dtype=int,
495 default=1,
496 )
497 ccdGrayFocalPlaneChebyshevOrder = pexConfig.Field(
498 doc="Order of the 2D chebyshev polynomials for focal plane fit.",
499 dtype=int,
500 default=3,
501 )
502 cycleNumber = pexConfig.Field(
503 doc=("FGCM fit cycle number. This is automatically incremented after each run "
504 "and stage of outlier rejection. See cookbook for details."),
505 dtype=int,
506 default=None,
507 )
508 isFinalCycle = pexConfig.Field(
509 doc=("Is this the final cycle of the fitting? Will automatically compute final "
510 "selection of stars and photometric exposures, and will output zeropoints "
511 "and standard stars for use in fgcmOutputProducts"),
512 dtype=bool,
513 default=False,
514 )
515 maxIterBeforeFinalCycle = pexConfig.Field(
516 doc=("Maximum fit iterations, prior to final cycle. The number of iterations "
517 "will always be 0 in the final cycle for cleanup and final selection."),
518 dtype=int,
519 default=50,
520 )
521 deltaMagBkgOffsetPercentile = pexConfig.Field(
522 doc=("Percentile brightest stars on a visit/ccd to use to compute net "
523 "offset from local background subtraction."),
524 dtype=float,
525 default=0.25,
526 )
527 deltaMagBkgPerCcd = pexConfig.Field(
528 doc=("Compute net offset from local background subtraction per-ccd? "
529 "Otherwise, use computation per visit."),
530 dtype=bool,
531 default=False,
532 )
533 utBoundary = pexConfig.Field(
534 doc="Boundary (in UTC) from day-to-day",
535 dtype=float,
536 default=None,
537 )
538 washMjds = pexConfig.ListField(
539 doc="Mirror wash MJDs",
540 dtype=float,
541 default=(0.0,),
542 )
543 epochMjds = pexConfig.ListField(
544 doc="Epoch boundaries in MJD",
545 dtype=float,
546 default=(0.0,),
547 )
548 minObsPerBand = pexConfig.Field(
549 doc="Minimum good observations per band",
550 dtype=int,
551 default=2,
552 )
553 # TODO: When DM-16511 is done, it will be possible to get the
554 # telescope latitude directly from the camera.
555 latitude = pexConfig.Field(
556 doc="Observatory latitude",
557 dtype=float,
558 default=None,
559 )
560 mirrorArea = pexConfig.Field(
561 doc="Mirror area (square meters) of telescope. If not set, will "
562 "try to estimate from camera.telescopeDiameter.",
563 dtype=float,
564 default=None,
565 optional=True,
566 )
567 defaultCameraOrientation = pexConfig.Field(
568 doc="Default camera orientation for QA plots.",
569 dtype=float,
570 default=None,
571 )
572 brightObsGrayMax = pexConfig.Field(
573 doc="Maximum gray extinction to be considered bright observation",
574 dtype=float,
575 default=0.15,
576 )
577 minStarPerCcd = pexConfig.Field(
578 doc=("Minimum number of good stars per CCD to be used in calibration fit. "
579 "CCDs with fewer stars will have their calibration estimated from other "
580 "CCDs in the same visit, with zeropoint error increased accordingly."),
581 dtype=int,
582 default=5,
583 )
584 minCcdPerExp = pexConfig.Field(
585 doc=("Minimum number of good CCDs per exposure/visit to be used in calibration fit. "
586 "Visits with fewer good CCDs will have CCD zeropoints estimated where possible."),
587 dtype=int,
588 default=5,
589 )
590 maxCcdGrayErr = pexConfig.Field(
591 doc="Maximum error on CCD gray offset to be considered photometric",
592 dtype=float,
593 default=0.05,
594 )
595 minStarPerExp = pexConfig.Field(
596 doc=("Minimum number of good stars per exposure/visit to be used in calibration fit. "
597 "Visits with fewer good stars will have CCD zeropoints estimated where possible."),
598 dtype=int,
599 default=600,
600 )
601 minExpPerNight = pexConfig.Field(
602 doc="Minimum number of good exposures/visits to consider a partly photometric night",
603 dtype=int,
604 default=10,
605 )
606 expGrayInitialCut = pexConfig.Field(
607 doc=("Maximum exposure/visit gray value for initial selection of possible photometric "
608 "observations."),
609 dtype=float,
610 default=-0.25,
611 )
612 expGrayPhotometricCutDict = pexConfig.DictField(
613 doc=("Per-band specification on maximum (negative) achromatic exposure residual "
614 "('gray term') for a visit to be considered photometric. Must have one "
615 "entry per band. Broad-band filters should be -0.05."),
616 keytype=str,
617 itemtype=float,
618 default={},
619 )
620 expGrayHighCutDict = pexConfig.DictField(
621 doc=("Per-band specification on maximum (positive) achromatic exposure residual "
622 "('gray term') for a visit to be considered photometric. Must have one "
623 "entry per band. Broad-band filters should be 0.2."),
624 keytype=str,
625 itemtype=float,
626 default={},
627 )
628 expGrayRecoverCut = pexConfig.Field(
629 doc=("Maximum (negative) exposure gray to be able to recover bad ccds via interpolation. "
630 "Visits with more gray extinction will only get CCD zeropoints if there are "
631 "sufficient star observations (minStarPerCcd) on that CCD."),
632 dtype=float,
633 default=-1.0,
634 )
635 expVarGrayPhotometricCutDict = pexConfig.DictField(
636 doc=("Per-band specification on maximum exposure variance to be considered possibly "
637 "photometric. Must have one entry per band. Broad-band filters should be "
638 "0.0005."),
639 keytype=str,
640 itemtype=float,
641 default={},
642 )
643 expGrayErrRecoverCut = pexConfig.Field(
644 doc=("Maximum exposure gray error to be able to recover bad ccds via interpolation. "
645 "Visits with more gray variance will only get CCD zeropoints if there are "
646 "sufficient star observations (minStarPerCcd) on that CCD."),
647 dtype=float,
648 default=0.05,
649 )
650 aperCorrFitNBins = pexConfig.Field(
651 doc=("Number of aperture bins used in aperture correction fit. When set to 0"
652 "no fit will be performed, and the config.aperCorrInputSlopes will be "
653 "used if available."),
654 dtype=int,
655 default=10,
656 )
657 aperCorrInputSlopeDict = pexConfig.DictField(
658 doc=("Per-band specification of aperture correction input slope parameters. These "
659 "are used on the first fit iteration, and aperture correction parameters will "
660 "be updated from the data if config.aperCorrFitNBins > 0. It is recommended "
661 "to set this when there is insufficient data to fit the parameters (e.g. "
662 "tract mode)."),
663 keytype=str,
664 itemtype=float,
665 default={},
666 )
667 sedboundaryterms = pexConfig.ConfigField(
668 doc="Mapping from bands to SED boundary term names used is sedterms.",
669 dtype=SedboundarytermDict,
670 )
671 sedterms = pexConfig.ConfigField(
672 doc="Mapping from terms to bands for fgcm linear SED approximations.",
673 dtype=SedtermDict,
674 )
675 sigFgcmMaxErr = pexConfig.Field(
676 doc="Maximum mag error for fitting sigma_FGCM",
677 dtype=float,
678 default=0.01,
679 )
680 sigFgcmMaxEGrayDict = pexConfig.DictField(
681 doc=("Per-band specification for maximum (absolute) achromatic residual (gray value) "
682 "for observations in sigma_fgcm (raw repeatability). Broad-band filters "
683 "should be 0.05."),
684 keytype=str,
685 itemtype=float,
686 default={},
687 )
688 ccdGrayMaxStarErr = pexConfig.Field(
689 doc=("Maximum error on a star observation to use in ccd gray (achromatic residual) "
690 "computation"),
691 dtype=float,
692 default=0.10,
693 )
694 approxThroughputDict = pexConfig.DictField(
695 doc=("Per-band specification of the approximate overall throughput at the start of "
696 "calibration observations. Must have one entry per band. Typically should "
697 "be 1.0."),
698 keytype=str,
699 itemtype=float,
700 default={},
701 )
702 sigmaCalRange = pexConfig.ListField(
703 doc="Allowed range for systematic error floor estimation",
704 dtype=float,
705 default=(0.001, 0.003),
706 )
707 sigmaCalFitPercentile = pexConfig.ListField(
708 doc="Magnitude percentile range to fit systematic error floor",
709 dtype=float,
710 default=(0.05, 0.15),
711 )
712 sigmaCalPlotPercentile = pexConfig.ListField(
713 doc="Magnitude percentile range to plot systematic error floor",
714 dtype=float,
715 default=(0.05, 0.95),
716 )
717 sigma0Phot = pexConfig.Field(
718 doc="Systematic error floor for all zeropoints",
719 dtype=float,
720 default=0.003,
721 )
722 mapLongitudeRef = pexConfig.Field(
723 doc="Reference longitude for plotting maps",
724 dtype=float,
725 default=0.0,
726 )
727 mapNSide = pexConfig.Field(
728 doc="Healpix nside for plotting maps",
729 dtype=int,
730 default=256,
731 )
732 outfileBase = pexConfig.Field(
733 doc="Filename start for plot output files",
734 dtype=str,
735 default=None,
736 )
737 starColorCuts = pexConfig.ListField(
738 doc=("Encoded star-color cuts (using calibration star colors). "
739 "This is a list with each entry a string of the format "
740 "``band1,band2,low,high`` such that only stars of color "
741 "low < band1 - band2 < high will be used for calibration."),
742 dtype=str,
743 default=("NO_DATA",),
744 )
745 refStarColorCuts = pexConfig.ListField(
746 doc=("Encoded star color cuts specifically to apply to reference stars. "
747 "This is a list with each entry a string of the format "
748 "``band1,band2,low,high`` such that only stars of color "
749 "low < band1 - band2 < high will be used as reference stars."),
750 dtype=str,
751 default=("NO_DATA",),
752 )
753 colorSplitBands = pexConfig.ListField(
754 doc="Band names to use to split stars by color. Must have 2 entries.",
755 dtype=str,
756 length=2,
757 default=('g', 'i'),
758 )
759 modelMagErrors = pexConfig.Field(
760 doc="Should FGCM model the magnitude errors from sky/fwhm? (False means trust inputs)",
761 dtype=bool,
762 default=True,
763 )
764 useQuadraticPwv = pexConfig.Field(
765 doc="Model PWV with a quadratic term for variation through the night?",
766 dtype=bool,
767 default=False,
768 )
769 instrumentParsPerBand = pexConfig.Field(
770 doc=("Model instrumental parameters per band? "
771 "Otherwise, instrumental parameters (QE changes with time) are "
772 "shared among all bands."),
773 dtype=bool,
774 default=False,
775 )
776 instrumentSlopeMinDeltaT = pexConfig.Field(
777 doc=("Minimum time change (in days) between observations to use in constraining "
778 "instrument slope."),
779 dtype=float,
780 default=20.0,
781 )
782 fitMirrorChromaticity = pexConfig.Field(
783 doc="Fit (intraband) mirror chromatic term?",
784 dtype=bool,
785 default=False,
786 )
787 fitCcdChromaticityDict = pexConfig.DictField(
788 doc="Specification on whether to compute first-order quantum efficiency (QE) "
789 "adjustments. Key is band, and value will be True or False. Any band "
790 "not explicitly specified will default to False.",
791 keytype=str,
792 itemtype=bool,
793 default={},
794 )
795 coatingMjds = pexConfig.ListField(
796 doc="Mirror coating dates in MJD",
797 dtype=float,
798 default=(0.0,),
799 )
800 outputStandardsBeforeFinalCycle = pexConfig.Field(
801 doc="Output standard stars prior to final cycle? Used in debugging.",
802 dtype=bool,
803 default=False,
804 )
805 outputZeropointsBeforeFinalCycle = pexConfig.Field(
806 doc="Output standard stars prior to final cycle? Used in debugging.",
807 dtype=bool,
808 default=False,
809 )
810 useRepeatabilityForExpGrayCutsDict = pexConfig.DictField(
811 doc=("Per-band specification on whether to use star repeatability (instead of exposures) "
812 "for computing photometric cuts. Recommended for tract mode or bands with few visits."),
813 keytype=str,
814 itemtype=bool,
815 default={},
816 )
817 autoPhotometricCutNSig = pexConfig.Field(
818 doc=("Number of sigma for automatic computation of (low) photometric cut. "
819 "Cut is based on exposure gray width (per band), unless "
820 "useRepeatabilityForExpGrayCuts is set, in which case the star "
821 "repeatability is used (also per band)."),
822 dtype=float,
823 default=3.0,
824 )
825 autoHighCutNSig = pexConfig.Field(
826 doc=("Number of sigma for automatic computation of (high) outlier cut. "
827 "Cut is based on exposure gray width (per band), unless "
828 "useRepeatabilityForExpGrayCuts is set, in which case the star "
829 "repeatability is used (also per band)."),
830 dtype=float,
831 default=4.0,
832 )
833 quietMode = pexConfig.Field(
834 doc="Be less verbose with logging.",
835 dtype=bool,
836 default=False,
837 )
838 doPlots = pexConfig.Field(
839 doc="Make fgcm QA plots.",
840 dtype=bool,
841 default=True,
842 )
843 randomSeed = pexConfig.Field(
844 doc="Random seed for fgcm for consistency in tests.",
845 dtype=int,
846 default=None,
847 optional=True,
848 )
849 deltaAperFitMinNgoodObs = pexConfig.Field(
850 doc="Minimum number of good observations to use mean delta-aper values in fits.",
851 dtype=int,
852 default=2,
853 )
854 deltaAperFitPerCcdNx = pexConfig.Field(
855 doc=("Number of x bins per ccd when computing delta-aper background offsets. "
856 "Only used when ``doComputeDeltaAperPerCcd`` is True."),
857 dtype=int,
858 default=10,
859 )
860 deltaAperFitPerCcdNy = pexConfig.Field(
861 doc=("Number of y bins per ccd when computing delta-aper background offsets. "
862 "Only used when ``doComputeDeltaAperPerCcd`` is True."),
863 dtype=int,
864 default=10,
865 )
866 deltaAperFitSpatialNside = pexConfig.Field(
867 doc="Healpix nside to compute spatial delta-aper background offset maps.",
868 dtype=int,
869 default=64,
870 )
871 deltaAperInnerRadiusArcsec = pexConfig.Field(
872 doc=("Inner radius used to compute deltaMagAper (arcseconds). "
873 "Must be positive and less than ``deltaAperOuterRadiusArcsec`` if "
874 "any of ``doComputeDeltaAperPerVisit``, ``doComputeDeltaAperPerStar``, "
875 "``doComputeDeltaAperMap``, ``doComputeDeltaAperPerCcd`` are set."),
876 dtype=float,
877 default=0.0,
878 )
879 deltaAperOuterRadiusArcsec = pexConfig.Field(
880 doc=("Outer radius used to compute deltaMagAper (arcseconds). "
881 "Must be positive and greater than ``deltaAperInnerRadiusArcsec`` if "
882 "any of ``doComputeDeltaAperPerVisit``, ``doComputeDeltaAperPerStar``, "
883 "``doComputeDeltaAperMap``, ``doComputeDeltaAperPerCcd`` are set."),
884 dtype=float,
885 default=0.0,
886 )
887 doComputeDeltaAperPerVisit = pexConfig.Field(
888 doc=("Do the computation of delta-aper background offsets per visit? "
889 "Note: this option can be very slow when there are many visits."),
890 dtype=bool,
891 default=False,
892 )
893 doComputeDeltaAperPerStar = pexConfig.Field(
894 doc="Do the computation of delta-aper mean values per star?",
895 dtype=bool,
896 default=True,
897 )
898 doComputeDeltaAperMap = pexConfig.Field(
899 doc=("Do the computation of delta-aper spatial maps? "
900 "This is only used if ``doComputeDeltaAperPerStar`` is True,"),
901 dtype=bool,
902 default=False,
903 )
904 doComputeDeltaAperPerCcd = pexConfig.Field(
905 doc="Do the computation of per-ccd delta-aper background offsets?",
906 dtype=bool,
907 default=False,
908 )
909
910 def validate(self):
911 super().validate()
912
913 if self.connections.previousCycleNumber != str(self.cycleNumber - 1):
914 msg = "cycleNumber in template must be connections.previousCycleNumber + 1"
915 raise RuntimeError(msg)
916 if self.connections.cycleNumber != str(self.cycleNumber):
917 msg = "cycleNumber in template must be equal to connections.cycleNumber"
918 raise RuntimeError(msg)
919
920 for band in self.fitBands:
921 if band not in self.bands:
922 msg = 'fitBand %s not in bands' % (band)
923 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.fitBands, self, msg)
924 for band in self.requiredBands:
925 if band not in self.bands:
926 msg = 'requiredBand %s not in bands' % (band)
927 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.requiredBands, self, msg)
928 for band in self.colorSplitBands:
929 if band not in self.bands:
930 msg = 'colorSplitBand %s not in bands' % (band)
931 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.colorSplitBands, self, msg)
932 for band in self.bands:
933 if band not in self.superStarSubCcdDict:
934 msg = 'band %s not in superStarSubCcdDict' % (band)
935 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.superStarSubCcdDict,
936 self, msg)
937 if band not in self.ccdGraySubCcdDict:
938 msg = 'band %s not in ccdGraySubCcdDict' % (band)
939 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.ccdGraySubCcdDict,
940 self, msg)
941 if band not in self.expGrayPhotometricCutDict:
942 msg = 'band %s not in expGrayPhotometricCutDict' % (band)
943 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayPhotometricCutDict,
944 self, msg)
945 if band not in self.expGrayHighCutDict:
946 msg = 'band %s not in expGrayHighCutDict' % (band)
947 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayHighCutDict,
948 self, msg)
949 if band not in self.expVarGrayPhotometricCutDict:
950 msg = 'band %s not in expVarGrayPhotometricCutDict' % (band)
951 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expVarGrayPhotometricCutDict,
952 self, msg)
953 if band not in self.sigFgcmMaxEGrayDict:
954 msg = 'band %s not in sigFgcmMaxEGrayDict' % (band)
955 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.sigFgcmMaxEGrayDict,
956 self, msg)
957 if band not in self.approxThroughputDict:
958 msg = 'band %s not in approxThroughputDict' % (band)
959 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.approxThroughputDict,
960 self, msg)
961 if band not in self.useRepeatabilityForExpGrayCutsDict:
962 msg = 'band %s not in useRepeatabilityForExpGrayCutsDict' % (band)
963 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict,
964 self, msg)
965
966 if self.doComputeDeltaAperPerVisit or self.doComputeDeltaAperMap \
967 or self.doComputeDeltaAperPerCcd:
968 if self.deltaAperInnerRadiusArcsec <= 0.0:
969 msg = 'deltaAperInnerRadiusArcsec must be positive if deltaAper computations are turned on.'
970 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperInnerRadiusArcsec,
971 self, msg)
972 if self.deltaAperOuterRadiusArcsec <= 0.0:
973 msg = 'deltaAperOuterRadiusArcsec must be positive if deltaAper computations are turned on.'
974 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperOuterRadiusArcsec,
975 self, msg)
976 if self.deltaAperOuterRadiusArcsec <= self.deltaAperInnerRadiusArcsec:
977 msg = ('deltaAperOuterRadiusArcsec must be greater than deltaAperInnerRadiusArcsec if '
978 'deltaAper computations are turned on.')
979 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperOuterRadiusArcsec,
980 self, msg)
981
982
983class FgcmFitCycleTask(pipeBase.PipelineTask):
984 """
985 Run Single fit cycle for FGCM global calibration
986 """
987
988 ConfigClass = FgcmFitCycleConfig
989 _DefaultName = "fgcmFitCycle"
990
991 def __init__(self, initInputs=None, **kwargs):
992 super().__init__(**kwargs)
993
994 def runQuantum(self, butlerQC, inputRefs, outputRefs):
995 camera = butlerQC.get(inputRefs.camera)
996
997 nCore = butlerQC.resources.num_cores
998
999 handleDict = {}
1000
1001 handleDict['fgcmLookUpTable'] = butlerQC.get(inputRefs.fgcmLookUpTable)
1002 handleDict['fgcmVisitCatalog'] = butlerQC.get(inputRefs.fgcmVisitCatalog)
1003
1004 if self.config.useParquetCatalogFormat:
1005 handleDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservationsParquet)
1006 handleDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIdsParquet)
1007 if self.config.doReferenceCalibration:
1008 handleDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStarsParquet)
1009 else:
1010 handleDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservations)
1011 handleDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIds)
1012 handleDict['fgcmStarIndices'] = butlerQC.get(inputRefs.fgcmStarIndices)
1013 if self.config.doReferenceCalibration:
1014 handleDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStars)
1015 if self.config.cycleNumber > 0:
1016 handleDict['fgcmFlaggedStars'] = butlerQC.get(inputRefs.fgcmFlaggedStarsInput)
1017 handleDict['fgcmFitParameters'] = butlerQC.get(inputRefs.fgcmFitParametersInput)
1018
1019 fgcmDatasetDict = None
1020 if self.config.doMultipleCycles:
1021 # Run multiple cycles at once.
1022 config = copy.copy(self.config)
1023 config.update(cycleNumber=0)
1024 for cycle in range(self.config.multipleCyclesFinalCycleNumber + 1):
1025 if cycle == self.config.multipleCyclesFinalCycleNumber:
1026 config.update(isFinalCycle=True)
1027
1028 if cycle > 0:
1029 handleDict['fgcmFlaggedStars'] = fgcmDatasetDict['fgcmFlaggedStars']
1030 handleDict['fgcmFitParameters'] = fgcmDatasetDict['fgcmFitParameters']
1031
1032 fgcmDatasetDict, config = self._fgcmFitCycle(camera, handleDict, config=config, nCore=nCore)
1033 butlerQC.put(fgcmDatasetDict['fgcmFitParameters'],
1034 getattr(outputRefs, f'fgcmFitParameters{cycle}'))
1035 butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'],
1036 getattr(outputRefs, f'fgcmFlaggedStars{cycle}'))
1037 if self.outputZeropoints:
1038 butlerQC.put(fgcmDatasetDict['fgcmZeropoints'],
1039 getattr(outputRefs, f'fgcmZeropoints{cycle}'))
1040 butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'],
1041 getattr(outputRefs, f'fgcmAtmosphereParameters{cycle}'))
1042 if self.outputStandards:
1043 butlerQC.put(fgcmDatasetDict['fgcmStandardStars'],
1044 getattr(outputRefs, f'fgcmStandardStars{cycle}'))
1045 else:
1046 # Run a single cycle
1047 fgcmDatasetDict, _ = self._fgcmFitCycle(camera, handleDict, nCore=nCore)
1048
1049 butlerQC.put(fgcmDatasetDict['fgcmFitParameters'], outputRefs.fgcmFitParameters)
1050 butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'], outputRefs.fgcmFlaggedStars)
1051 if self.outputZeropoints:
1052 butlerQC.put(fgcmDatasetDict['fgcmZeropoints'], outputRefs.fgcmZeropoints)
1053 butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'], outputRefs.fgcmAtmosphereParameters)
1054 if self.outputStandards:
1055 butlerQC.put(fgcmDatasetDict['fgcmStandardStars'], outputRefs.fgcmStandardStars)
1056
1057 def _fgcmFitCycle(self, camera, handleDict, config=None, nCore=1):
1058 """
1059 Run the fit cycle
1060
1061 Parameters
1062 ----------
1063 camera : `lsst.afw.cameraGeom.Camera`
1064 handleDict : `dict`
1065 All handles are `lsst.daf.butler.DeferredDatasetHandle`
1066 handle dictionary with keys:
1067
1068 ``"fgcmLookUpTable"``
1069 handle for the FGCM look-up table.
1070 ``"fgcmVisitCatalog"``
1071 handle for visit summary catalog.
1072 ``"fgcmStarObservations"``
1073 handle for star observation catalog.
1074 ``"fgcmStarIds"``
1075 handle for star id catalog.
1076 ``"fgcmStarIndices"``
1077 handle for star index catalog.
1078 ``"fgcmReferenceStars"``
1079 handle for matched reference star catalog.
1080 ``"fgcmFlaggedStars"``
1081 handle for flagged star catalog.
1082 ``"fgcmFitParameters"``
1083 handle for fit parameter catalog.
1084 config : `lsst.pex.config.Config`, optional
1085 Configuration to use to override self.config.
1086 nCore : `int`, optional
1087 Number of cores to use during fitting.
1088
1089 Returns
1090 -------
1091 fgcmDatasetDict : `dict`
1092 Dictionary of datasets to persist.
1093 """
1094 if config is not None:
1095 _config = config
1096 else:
1097 _config = self.config
1098
1099 # Set defaults on whether to output standards and zeropoints
1100 self.maxIter = _config.maxIterBeforeFinalCycle
1101 self.outputStandards = _config.outputStandardsBeforeFinalCycle
1102 self.outputZeropoints = _config.outputZeropointsBeforeFinalCycle
1103 self.resetFitParameters = True
1104
1105 if _config.isFinalCycle:
1106 # This is the final fit cycle, so we do not want to reset fit
1107 # parameters, we want to run a final "clean-up" with 0 fit iterations,
1108 # and we always want to output standards and zeropoints
1109 self.maxIter = 0
1110 self.outputStandards = True
1111 self.outputZeropoints = True
1112 self.resetFitParameters = False
1113
1114 lutCat = handleDict['fgcmLookUpTable'].get()
1115 fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat,
1116 dict(_config.physicalFilterMap))
1117 del lutCat
1118
1119 configDict = makeConfigDict(_config, self.log, camera,
1120 self.maxIter, self.resetFitParameters,
1121 self.outputZeropoints,
1122 lutIndexVals[0]['FILTERNAMES'],
1123 nCore=nCore)
1124
1125 # next we need the exposure/visit information
1126 visitCat = handleDict['fgcmVisitCatalog'].get()
1127 fgcmExpInfo = translateVisitCatalog(visitCat)
1128 del visitCat
1129
1130 focalPlaneProjector = FocalPlaneProjector(camera,
1131 self.config.defaultCameraOrientation)
1132
1133 noFitsDict = {'lutIndex': lutIndexVals,
1134 'lutStd': lutStd,
1135 'expInfo': fgcmExpInfo,
1136 'focalPlaneProjector': focalPlaneProjector}
1137
1138 # set up the fitter object
1139 fgcmFitCycle = fgcm.FgcmFitCycle(configDict, useFits=False,
1140 noFitsDict=noFitsDict, noOutput=True)
1141
1142 # create the parameter object
1143 if (fgcmFitCycle.initialCycle):
1144 # cycle = 0, initial cycle
1145 fgcmPars = fgcm.FgcmParameters.newParsWithArrays(fgcmFitCycle.fgcmConfig,
1146 fgcmLut,
1147 fgcmExpInfo)
1148 else:
1149 if isinstance(handleDict['fgcmFitParameters'], afwTable.BaseCatalog):
1150 parCat = handleDict['fgcmFitParameters']
1151 else:
1152 parCat = handleDict['fgcmFitParameters'].get()
1153 inParInfo, inParams, inSuperStar = self._loadParameters(parCat)
1154 del parCat
1155 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
1156 fgcmExpInfo,
1157 inParInfo,
1158 inParams,
1159 inSuperStar)
1160
1161 # set up the stars...
1162 fgcmStars = fgcm.FgcmStars(fgcmFitCycle.fgcmConfig)
1163
1164 starObs = handleDict['fgcmStarObservations'].get()
1165 starIds = handleDict['fgcmStarIds'].get()
1166 if not self.config.useParquetCatalogFormat:
1167 starIndices = handleDict['fgcmStarIndices'].get()
1168 else:
1169 starIndices = None
1170
1171 # grab the flagged stars if available
1172 if 'fgcmFlaggedStars' in handleDict:
1173 if isinstance(handleDict['fgcmFlaggedStars'], afwTable.BaseCatalog):
1174 flaggedStars = handleDict['fgcmFlaggedStars']
1175 else:
1176 flaggedStars = handleDict['fgcmFlaggedStars'].get()
1177 flagId = flaggedStars['objId'][:]
1178 flagFlag = flaggedStars['objFlag'][:]
1179
1180 del flaggedStars
1181 elif self.config.useParquetCatalogFormat:
1182 # If we are using the parquet catalog format, then that means that
1183 # reserved stars have already been flagged. We extract the flags here
1184 # to input to fgcm, which will then be persisted (with additional
1185 # quality flags) as the fgcmFlaggedStars datatype in subsequent
1186 # fit cycles.
1187 (flagged,) = (starIds['obj_flag'] > 0).nonzero()
1188 flagId = starIds['fgcm_id'][flagged]
1189 flagFlag = starIds['obj_flag'][flagged]
1190 else:
1191 flagId = None
1192 flagFlag = None
1193
1194 if _config.doReferenceCalibration:
1195 refStars = handleDict['fgcmReferenceStars'].get()
1196
1197 refMag, refMagErr = extractReferenceMags(refStars,
1198 _config.bands,
1199 _config.physicalFilterMap)
1200
1201 refId = refStars['fgcm_id'][:]
1202 else:
1203 refStars = None
1204 refId = None
1205 refMag = None
1206 refMagErr = None
1207
1208 # match star observations to visits
1209 # Only those star observations that match visits from fgcmExpInfo['VISIT'] will
1210 # actually be transferred into fgcm using the indexing below.
1211 if self.config.useParquetCatalogFormat:
1212 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit'])
1213 else:
1214 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit'][starIndices['obsIndex']])
1215
1216 # The fgcmStars.loadStars method will copy all the star information into
1217 # special shared memory objects that will not blow up the memory usage when
1218 # used with python multiprocessing. Once all the numbers are copied,
1219 # it is necessary to release all references to the objects that previously
1220 # stored the data to ensure that the garbage collector can clear the memory,
1221 # and ensure that this memory is not copied when multiprocessing kicks in.
1222
1223 if self.config.useParquetCatalogFormat:
1224 # Note that the ra/dec coordinates for the parquet format are in
1225 # degrees, which is what fgcm expects.
1226 fgcmStars.loadStars(fgcmPars,
1227 starObs['visit'][:],
1228 starObs['detector'][:],
1229 starObs['ra'][:],
1230 starObs['dec'][:],
1231 starObs['inst_mag'][:],
1232 starObs['inst_mag_err'][:],
1233 fgcmExpInfo['FILTERNAME'][visitIndex],
1234 starIds['fgcm_id'][:],
1235 starIds['ra'][:],
1236 starIds['dec'][:],
1237 starIds['obs_arr_index'][:],
1238 starIds['n_obs'][:],
1239 obsX=starObs['x'][:],
1240 obsY=starObs['y'][:],
1241 obsDeltaMagBkg=starObs['delta_mag_bkg'][:],
1242 obsDeltaAper=starObs['delta_mag_aper'][:],
1243 refID=refId,
1244 refMag=refMag,
1245 refMagErr=refMagErr,
1246 flagID=flagId,
1247 flagFlag=flagFlag,
1248 computeNobs=True)
1249 else:
1250 # We determine the conversion from the native units (typically radians) to
1251 # degrees for the first star. This allows us to treat coord_ra/coord_dec as
1252 # numpy arrays rather than Angles, which would we approximately 600x slower.
1253 conv = starObs[0]['ra'].asDegrees() / float(starObs[0]['ra'])
1254
1255 fgcmStars.loadStars(fgcmPars,
1256 starObs['visit'][starIndices['obsIndex']],
1257 starObs['ccd'][starIndices['obsIndex']],
1258 starObs['ra'][starIndices['obsIndex']] * conv,
1259 starObs['dec'][starIndices['obsIndex']] * conv,
1260 starObs['instMag'][starIndices['obsIndex']],
1261 starObs['instMagErr'][starIndices['obsIndex']],
1262 fgcmExpInfo['FILTERNAME'][visitIndex],
1263 starIds['fgcm_id'][:],
1264 starIds['ra'][:],
1265 starIds['dec'][:],
1266 starIds['obsArrIndex'][:],
1267 starIds['nObs'][:],
1268 obsX=starObs['x'][starIndices['obsIndex']],
1269 obsY=starObs['y'][starIndices['obsIndex']],
1270 obsDeltaMagBkg=starObs['deltaMagBkg'][starIndices['obsIndex']],
1271 obsDeltaAper=starObs['deltaMagAper'][starIndices['obsIndex']],
1272 psfCandidate=starObs['psf_candidate'][starIndices['obsIndex']],
1273 refID=refId,
1274 refMag=refMag,
1275 refMagErr=refMagErr,
1276 flagID=flagId,
1277 flagFlag=flagFlag,
1278 computeNobs=True)
1279
1280 # Release all references to temporary objects holding star data (see above)
1281 del starObs
1282 del starIds
1283 del starIndices
1284 del flagId
1285 del flagFlag
1286 del refStars
1287 del refId
1288 del refMag
1289 del refMagErr
1290
1291 # and set the bits in the cycle object
1292 fgcmFitCycle.setLUT(fgcmLut)
1293 fgcmFitCycle.setStars(fgcmStars, fgcmPars)
1294 fgcmFitCycle.setPars(fgcmPars)
1295
1296 # finish the setup
1297 fgcmFitCycle.finishSetup()
1298
1299 # and run
1300 fgcmFitCycle.run()
1301
1302
1305
1306 fgcmDatasetDict = self._makeFgcmOutputDatasets(fgcmFitCycle)
1307
1308 # Output the config for the next cycle
1309 # We need to make a copy since the input one has been frozen
1310
1311 updatedPhotometricCutDict = {b: float(fgcmFitCycle.updatedPhotometricCut[i]) for
1312 i, b in enumerate(_config.bands)}
1313 updatedHighCutDict = {band: float(fgcmFitCycle.updatedHighCut[i]) for
1314 i, band in enumerate(_config.bands)}
1315
1316 outConfig = copy.copy(_config)
1317 outConfig.update(cycleNumber=(_config.cycleNumber + 1),
1318 precomputeSuperStarInitialCycle=False,
1319 freezeStdAtmosphere=False,
1320 expGrayPhotometricCutDict=updatedPhotometricCutDict,
1321 expGrayHighCutDict=updatedHighCutDict)
1322
1323 outConfig.connections.update(previousCycleNumber=str(_config.cycleNumber),
1324 cycleNumber=str(_config.cycleNumber + 1))
1325
1326 configFileName = '%s_cycle%02d_config.py' % (outConfig.outfileBase,
1327 outConfig.cycleNumber)
1328 outConfig.save(configFileName)
1329
1330 if _config.isFinalCycle == 1:
1331 # We are done, ready to output products
1332 self.log.info("Everything is in place to run fgcmOutputProducts.py")
1333 else:
1334 self.log.info("Saved config for next cycle to %s" % (configFileName))
1335 self.log.info("Be sure to look at:")
1336 self.log.info(" config.expGrayPhotometricCut")
1337 self.log.info(" config.expGrayHighCut")
1338 self.log.info("If you are satisfied with the fit, please set:")
1339 self.log.info(" config.isFinalCycle = True")
1340
1341 fgcmFitCycle.freeSharedMemory()
1342
1343 return fgcmDatasetDict, outConfig
1344
1345 def _loadParameters(self, parCat):
1346 """
1347 Load FGCM parameters from a previous fit cycle
1348
1349 Parameters
1350 ----------
1351 parCat : `lsst.afw.table.BaseCatalog`
1352 Parameter catalog in afw table form.
1353
1354 Returns
1355 -------
1356 inParInfo: `numpy.ndarray`
1357 Numpy array parameter information formatted for input to fgcm
1358 inParameters: `numpy.ndarray`
1359 Numpy array parameter values formatted for input to fgcm
1360 inSuperStar: `numpy.array`
1361 Superstar flat formatted for input to fgcm
1362 """
1363 parLutFilterNames = np.array(parCat[0]['lutFilterNames'].split(','))
1364 parFitBands = np.array(parCat[0]['fitBands'].split(','))
1365
1366 inParInfo = np.zeros(1, dtype=[('NCCD', 'i4'),
1367 ('LUTFILTERNAMES', parLutFilterNames.dtype.str,
1368 (parLutFilterNames.size, )),
1369 ('FITBANDS', parFitBands.dtype.str, (parFitBands.size, )),
1370 ('LNTAUUNIT', 'f8'),
1371 ('LNTAUSLOPEUNIT', 'f8'),
1372 ('ALPHAUNIT', 'f8'),
1373 ('LNPWVUNIT', 'f8'),
1374 ('LNPWVSLOPEUNIT', 'f8'),
1375 ('LNPWVQUADRATICUNIT', 'f8'),
1376 ('LNPWVGLOBALUNIT', 'f8'),
1377 ('O3UNIT', 'f8'),
1378 ('QESYSUNIT', 'f8'),
1379 ('FILTEROFFSETUNIT', 'f8'),
1380 ('HASEXTERNALPWV', 'i2'),
1381 ('HASEXTERNALTAU', 'i2')])
1382 inParInfo['NCCD'] = parCat['nCcd']
1383 inParInfo['LUTFILTERNAMES'][:] = parLutFilterNames
1384 inParInfo['FITBANDS'][:] = parFitBands
1385 inParInfo['HASEXTERNALPWV'] = parCat['hasExternalPwv']
1386 inParInfo['HASEXTERNALTAU'] = parCat['hasExternalTau']
1387
1388 inParams = np.zeros(1, dtype=[('PARALPHA', 'f8', (parCat['parAlpha'].size, )),
1389 ('PARO3', 'f8', (parCat['parO3'].size, )),
1390 ('PARLNTAUINTERCEPT', 'f8',
1391 (parCat['parLnTauIntercept'].size, )),
1392 ('PARLNTAUSLOPE', 'f8',
1393 (parCat['parLnTauSlope'].size, )),
1394 ('PARLNPWVINTERCEPT', 'f8',
1395 (parCat['parLnPwvIntercept'].size, )),
1396 ('PARLNPWVSLOPE', 'f8',
1397 (parCat['parLnPwvSlope'].size, )),
1398 ('PARLNPWVQUADRATIC', 'f8',
1399 (parCat['parLnPwvQuadratic'].size, )),
1400 ('PARQESYSINTERCEPT', 'f8',
1401 (parCat['parQeSysIntercept'].size, )),
1402 ('COMPQESYSSLOPE', 'f8',
1403 (parCat['compQeSysSlope'].size, )),
1404 ('PARFILTEROFFSET', 'f8',
1405 (parCat['parFilterOffset'].size, )),
1406 ('PARFILTEROFFSETFITFLAG', 'i2',
1407 (parCat['parFilterOffsetFitFlag'].size, )),
1408 ('PARRETRIEVEDLNPWVSCALE', 'f8'),
1409 ('PARRETRIEVEDLNPWVOFFSET', 'f8'),
1410 ('PARRETRIEVEDLNPWVNIGHTLYOFFSET', 'f8',
1411 (parCat['parRetrievedLnPwvNightlyOffset'].size, )),
1412 ('COMPABSTHROUGHPUT', 'f8',
1413 (parCat['compAbsThroughput'].size, )),
1414 ('COMPREFOFFSET', 'f8',
1415 (parCat['compRefOffset'].size, )),
1416 ('COMPREFSIGMA', 'f8',
1417 (parCat['compRefSigma'].size, )),
1418 ('COMPMIRRORCHROMATICITY', 'f8',
1419 (parCat['compMirrorChromaticity'].size, )),
1420 ('MIRRORCHROMATICITYPIVOT', 'f8',
1421 (parCat['mirrorChromaticityPivot'].size, )),
1422 ('COMPCCDCHROMATICITY', 'f8',
1423 (parCat['compCcdChromaticity'].size, )),
1424 ('COMPMEDIANSEDSLOPE', 'f8',
1425 (parCat['compMedianSedSlope'].size, )),
1426 ('COMPAPERCORRPIVOT', 'f8',
1427 (parCat['compAperCorrPivot'].size, )),
1428 ('COMPAPERCORRSLOPE', 'f8',
1429 (parCat['compAperCorrSlope'].size, )),
1430 ('COMPAPERCORRSLOPEERR', 'f8',
1431 (parCat['compAperCorrSlopeErr'].size, )),
1432 ('COMPAPERCORRRANGE', 'f8',
1433 (parCat['compAperCorrRange'].size, )),
1434 ('COMPMODELERREXPTIMEPIVOT', 'f8',
1435 (parCat['compModelErrExptimePivot'].size, )),
1436 ('COMPMODELERRFWHMPIVOT', 'f8',
1437 (parCat['compModelErrFwhmPivot'].size, )),
1438 ('COMPMODELERRSKYPIVOT', 'f8',
1439 (parCat['compModelErrSkyPivot'].size, )),
1440 ('COMPMODELERRPARS', 'f8',
1441 (parCat['compModelErrPars'].size, )),
1442 ('COMPEXPGRAY', 'f8',
1443 (parCat['compExpGray'].size, )),
1444 ('COMPVARGRAY', 'f8',
1445 (parCat['compVarGray'].size, )),
1446 ('COMPEXPDELTAMAGBKG', 'f8',
1447 (parCat['compExpDeltaMagBkg'].size, )),
1448 ('COMPNGOODSTARPEREXP', 'i4',
1449 (parCat['compNGoodStarPerExp'].size, )),
1450 ('COMPEXPREFOFFSET', 'f8',
1451 (parCat['compExpRefOffset'].size, )),
1452 ('COMPSIGFGCM', 'f8',
1453 (parCat['compSigFgcm'].size, )),
1454 ('COMPSIGMACAL', 'f8',
1455 (parCat['compSigmaCal'].size, )),
1456 ('COMPRETRIEVEDLNPWV', 'f8',
1457 (parCat['compRetrievedLnPwv'].size, )),
1458 ('COMPRETRIEVEDLNPWVRAW', 'f8',
1459 (parCat['compRetrievedLnPwvRaw'].size, )),
1460 ('COMPRETRIEVEDLNPWVFLAG', 'i2',
1461 (parCat['compRetrievedLnPwvFlag'].size, )),
1462 ('COMPRETRIEVEDTAUNIGHT', 'f8',
1463 (parCat['compRetrievedTauNight'].size, )),
1464 ('COMPEPSILON', 'f8',
1465 (parCat['compEpsilon'].size, )),
1466 ('COMPMEDDELTAAPER', 'f8',
1467 (parCat['compMedDeltaAper'].size, )),
1468 ('COMPGLOBALEPSILON', 'f4',
1469 (parCat['compGlobalEpsilon'].size, )),
1470 ('COMPEPSILONMAP', 'f4',
1471 (parCat['compEpsilonMap'].size, )),
1472 ('COMPEPSILONNSTARMAP', 'i4',
1473 (parCat['compEpsilonNStarMap'].size, )),
1474 ('COMPEPSILONCCDMAP', 'f4',
1475 (parCat['compEpsilonCcdMap'].size, )),
1476 ('COMPEPSILONCCDNSTARMAP', 'i4',
1477 (parCat['compEpsilonCcdNStarMap'].size, ))])
1478
1479 inParams['PARALPHA'][:] = parCat['parAlpha'][0, :]
1480 inParams['PARO3'][:] = parCat['parO3'][0, :]
1481 inParams['PARLNTAUINTERCEPT'][:] = parCat['parLnTauIntercept'][0, :]
1482 inParams['PARLNTAUSLOPE'][:] = parCat['parLnTauSlope'][0, :]
1483 inParams['PARLNPWVINTERCEPT'][:] = parCat['parLnPwvIntercept'][0, :]
1484 inParams['PARLNPWVSLOPE'][:] = parCat['parLnPwvSlope'][0, :]
1485 inParams['PARLNPWVQUADRATIC'][:] = parCat['parLnPwvQuadratic'][0, :]
1486 inParams['PARQESYSINTERCEPT'][:] = parCat['parQeSysIntercept'][0, :]
1487 inParams['COMPQESYSSLOPE'][:] = parCat['compQeSysSlope'][0, :]
1488 inParams['PARFILTEROFFSET'][:] = parCat['parFilterOffset'][0, :]
1489 inParams['PARFILTEROFFSETFITFLAG'][:] = parCat['parFilterOffsetFitFlag'][0, :]
1490 inParams['PARRETRIEVEDLNPWVSCALE'] = parCat['parRetrievedLnPwvScale']
1491 inParams['PARRETRIEVEDLNPWVOFFSET'] = parCat['parRetrievedLnPwvOffset']
1492 inParams['PARRETRIEVEDLNPWVNIGHTLYOFFSET'][:] = parCat['parRetrievedLnPwvNightlyOffset'][0, :]
1493 inParams['COMPABSTHROUGHPUT'][:] = parCat['compAbsThroughput'][0, :]
1494 inParams['COMPREFOFFSET'][:] = parCat['compRefOffset'][0, :]
1495 inParams['COMPREFSIGMA'][:] = parCat['compRefSigma'][0, :]
1496 inParams['COMPMIRRORCHROMATICITY'][:] = parCat['compMirrorChromaticity'][0, :]
1497 inParams['MIRRORCHROMATICITYPIVOT'][:] = parCat['mirrorChromaticityPivot'][0, :]
1498 inParams['COMPCCDCHROMATICITY'][:] = parCat['compCcdChromaticity'][0, :]
1499 inParams['COMPMEDIANSEDSLOPE'][:] = parCat['compMedianSedSlope'][0, :]
1500 inParams['COMPAPERCORRPIVOT'][:] = parCat['compAperCorrPivot'][0, :]
1501 inParams['COMPAPERCORRSLOPE'][:] = parCat['compAperCorrSlope'][0, :]
1502 inParams['COMPAPERCORRSLOPEERR'][:] = parCat['compAperCorrSlopeErr'][0, :]
1503 inParams['COMPAPERCORRRANGE'][:] = parCat['compAperCorrRange'][0, :]
1504 inParams['COMPMODELERREXPTIMEPIVOT'][:] = parCat['compModelErrExptimePivot'][0, :]
1505 inParams['COMPMODELERRFWHMPIVOT'][:] = parCat['compModelErrFwhmPivot'][0, :]
1506 inParams['COMPMODELERRSKYPIVOT'][:] = parCat['compModelErrSkyPivot'][0, :]
1507 inParams['COMPMODELERRPARS'][:] = parCat['compModelErrPars'][0, :]
1508 inParams['COMPEXPGRAY'][:] = parCat['compExpGray'][0, :]
1509 inParams['COMPVARGRAY'][:] = parCat['compVarGray'][0, :]
1510 inParams['COMPEXPDELTAMAGBKG'][:] = parCat['compExpDeltaMagBkg'][0, :]
1511 inParams['COMPNGOODSTARPEREXP'][:] = parCat['compNGoodStarPerExp'][0, :]
1512 inParams['COMPEXPREFOFFSET'][:] = parCat['compExpRefOffset'][0, :]
1513 inParams['COMPSIGFGCM'][:] = parCat['compSigFgcm'][0, :]
1514 inParams['COMPSIGMACAL'][:] = parCat['compSigmaCal'][0, :]
1515 inParams['COMPRETRIEVEDLNPWV'][:] = parCat['compRetrievedLnPwv'][0, :]
1516 inParams['COMPRETRIEVEDLNPWVRAW'][:] = parCat['compRetrievedLnPwvRaw'][0, :]
1517 inParams['COMPRETRIEVEDLNPWVFLAG'][:] = parCat['compRetrievedLnPwvFlag'][0, :]
1518 inParams['COMPRETRIEVEDTAUNIGHT'][:] = parCat['compRetrievedTauNight'][0, :]
1519 inParams['COMPEPSILON'][:] = parCat['compEpsilon'][0, :]
1520 inParams['COMPMEDDELTAAPER'][:] = parCat['compMedDeltaAper'][0, :]
1521 inParams['COMPGLOBALEPSILON'][:] = parCat['compGlobalEpsilon'][0, :]
1522 inParams['COMPEPSILONMAP'][:] = parCat['compEpsilonMap'][0, :]
1523 inParams['COMPEPSILONNSTARMAP'][:] = parCat['compEpsilonNStarMap'][0, :]
1524 inParams['COMPEPSILONCCDMAP'][:] = parCat['compEpsilonCcdMap'][0, :]
1525 inParams['COMPEPSILONCCDNSTARMAP'][:] = parCat['compEpsilonCcdNStarMap'][0, :]
1526
1527 inSuperStar = np.zeros(parCat['superstarSize'][0, :], dtype='f8')
1528 inSuperStar[:, :, :, :] = parCat['superstar'][0, :].reshape(inSuperStar.shape)
1529
1530 return (inParInfo, inParams, inSuperStar)
1531
1532 def _makeFgcmOutputDatasets(self, fgcmFitCycle):
1533 """
1534 Persist FGCM datasets through the butler.
1535
1536 Parameters
1537 ----------
1538 fgcmFitCycle: `lsst.fgcm.FgcmFitCycle`
1539 Fgcm Fit cycle object
1540 """
1541 fgcmDatasetDict = {}
1542
1543 # Save the parameters
1544 parInfo, pars = fgcmFitCycle.fgcmPars.parsToArrays()
1545
1546 parSchema = afwTable.Schema()
1547
1548 comma = ','
1549 lutFilterNameString = comma.join([n.decode('utf-8')
1550 for n in parInfo['LUTFILTERNAMES'][0]])
1551 fitBandString = comma.join([n.decode('utf-8')
1552 for n in parInfo['FITBANDS'][0]])
1553
1554 parSchema = self._makeParSchema(parInfo, pars, fgcmFitCycle.fgcmPars.parSuperStarFlat,
1555 lutFilterNameString, fitBandString)
1556 parCat = self._makeParCatalog(parSchema, parInfo, pars,
1557 fgcmFitCycle.fgcmPars.parSuperStarFlat,
1558 lutFilterNameString, fitBandString)
1559
1560 fgcmDatasetDict['fgcmFitParameters'] = parCat
1561
1562 # Save the indices of the flagged stars
1563 # (stars that have been (a) reserved from the fit for testing and
1564 # (b) bad stars that have failed quality checks.)
1565 flagStarSchema = self._makeFlagStarSchema()
1566 flagStarStruct = fgcmFitCycle.fgcmStars.getFlagStarIndices()
1567 flagStarCat = self._makeFlagStarCat(flagStarSchema, flagStarStruct)
1568
1569 fgcmDatasetDict['fgcmFlaggedStars'] = flagStarCat
1570
1571 # Save the zeropoint information and atmospheres only if desired
1572 if self.outputZeropoints:
1573 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
1574 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
1575
1576 zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
1577 zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
1578
1579 fgcmDatasetDict['fgcmZeropoints'] = zptCat
1580
1581 # Save atmosphere values
1582 # These are generated by the same code that generates zeropoints
1583 atmSchema = makeAtmSchema()
1584 atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
1585
1586 fgcmDatasetDict['fgcmAtmosphereParameters'] = atmCat
1587
1588 # Save the standard stars (if configured)
1589 if self.outputStandards:
1590 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
1591 stdSchema = makeStdSchema(len(goodBands))
1592 stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
1593
1594 fgcmDatasetDict['fgcmStandardStars'] = stdCat
1595
1596 return fgcmDatasetDict
1597
1598 def _makeParSchema(self, parInfo, pars, parSuperStarFlat,
1599 lutFilterNameString, fitBandString):
1600 """
1601 Make the parameter persistence schema
1602
1603 Parameters
1604 ----------
1605 parInfo: `numpy.ndarray`
1606 Parameter information returned by fgcm
1607 pars: `numpy.ndarray`
1608 Parameter values returned by fgcm
1609 parSuperStarFlat: `numpy.array`
1610 Superstar flat values returned by fgcm
1611 lutFilterNameString: `str`
1612 Combined string of all the lutFilterNames
1613 fitBandString: `str`
1614 Combined string of all the fitBands
1615
1616 Returns
1617 -------
1618 parSchema: `afwTable.schema`
1619 """
1620
1621 parSchema = afwTable.Schema()
1622
1623 # parameter info section
1624 parSchema.addField('nCcd', type=np.int32, doc='Number of CCDs')
1625 parSchema.addField('lutFilterNames', type=str, doc='LUT Filter names in parameter file',
1626 size=len(lutFilterNameString))
1627 parSchema.addField('fitBands', type=str, doc='Bands that were fit',
1628 size=len(fitBandString))
1629 parSchema.addField('lnTauUnit', type=np.float64, doc='Step units for ln(AOD)')
1630 parSchema.addField('lnTauSlopeUnit', type=np.float64,
1631 doc='Step units for ln(AOD) slope')
1632 parSchema.addField('alphaUnit', type=np.float64, doc='Step units for alpha')
1633 parSchema.addField('lnPwvUnit', type=np.float64, doc='Step units for ln(pwv)')
1634 parSchema.addField('lnPwvSlopeUnit', type=np.float64,
1635 doc='Step units for ln(pwv) slope')
1636 parSchema.addField('lnPwvQuadraticUnit', type=np.float64,
1637 doc='Step units for ln(pwv) quadratic term')
1638 parSchema.addField('lnPwvGlobalUnit', type=np.float64,
1639 doc='Step units for global ln(pwv) parameters')
1640 parSchema.addField('o3Unit', type=np.float64, doc='Step units for O3')
1641 parSchema.addField('qeSysUnit', type=np.float64, doc='Step units for mirror gray')
1642 parSchema.addField('filterOffsetUnit', type=np.float64, doc='Step units for filter offset')
1643 parSchema.addField('hasExternalPwv', type=np.int32, doc='Parameters fit using external pwv')
1644 parSchema.addField('hasExternalTau', type=np.int32, doc='Parameters fit using external tau')
1645
1646 # parameter section
1647 parSchema.addField('parAlpha', type='ArrayD', doc='Alpha parameter vector',
1648 size=pars['PARALPHA'].size)
1649 parSchema.addField('parO3', type='ArrayD', doc='O3 parameter vector',
1650 size=pars['PARO3'].size)
1651 parSchema.addField('parLnTauIntercept', type='ArrayD',
1652 doc='ln(Tau) intercept parameter vector',
1653 size=pars['PARLNTAUINTERCEPT'].size)
1654 parSchema.addField('parLnTauSlope', type='ArrayD',
1655 doc='ln(Tau) slope parameter vector',
1656 size=pars['PARLNTAUSLOPE'].size)
1657 parSchema.addField('parLnPwvIntercept', type='ArrayD', doc='ln(pwv) intercept parameter vector',
1658 size=pars['PARLNPWVINTERCEPT'].size)
1659 parSchema.addField('parLnPwvSlope', type='ArrayD', doc='ln(pwv) slope parameter vector',
1660 size=pars['PARLNPWVSLOPE'].size)
1661 parSchema.addField('parLnPwvQuadratic', type='ArrayD', doc='ln(pwv) quadratic parameter vector',
1662 size=pars['PARLNPWVQUADRATIC'].size)
1663 parSchema.addField('parQeSysIntercept', type='ArrayD', doc='Mirror gray intercept parameter vector',
1664 size=pars['PARQESYSINTERCEPT'].size)
1665 parSchema.addField('compQeSysSlope', type='ArrayD', doc='Mirror gray slope parameter vector',
1666 size=pars[0]['COMPQESYSSLOPE'].size)
1667 parSchema.addField('parFilterOffset', type='ArrayD', doc='Filter offset parameter vector',
1668 size=pars['PARFILTEROFFSET'].size)
1669 parSchema.addField('parFilterOffsetFitFlag', type='ArrayI', doc='Filter offset parameter fit flag',
1670 size=pars['PARFILTEROFFSETFITFLAG'].size)
1671 parSchema.addField('parRetrievedLnPwvScale', type=np.float64,
1672 doc='Global scale for retrieved ln(pwv)')
1673 parSchema.addField('parRetrievedLnPwvOffset', type=np.float64,
1674 doc='Global offset for retrieved ln(pwv)')
1675 parSchema.addField('parRetrievedLnPwvNightlyOffset', type='ArrayD',
1676 doc='Nightly offset for retrieved ln(pwv)',
1677 size=pars['PARRETRIEVEDLNPWVNIGHTLYOFFSET'].size)
1678 parSchema.addField('compAbsThroughput', type='ArrayD',
1679 doc='Absolute throughput (relative to transmission curves)',
1680 size=pars['COMPABSTHROUGHPUT'].size)
1681 parSchema.addField('compRefOffset', type='ArrayD',
1682 doc='Offset between reference stars and calibrated stars',
1683 size=pars['COMPREFOFFSET'].size)
1684 parSchema.addField('compRefSigma', type='ArrayD',
1685 doc='Width of reference star/calibrated star distribution',
1686 size=pars['COMPREFSIGMA'].size)
1687 parSchema.addField('compMirrorChromaticity', type='ArrayD',
1688 doc='Computed mirror chromaticity terms',
1689 size=pars['COMPMIRRORCHROMATICITY'].size)
1690 parSchema.addField('mirrorChromaticityPivot', type='ArrayD',
1691 doc='Mirror chromaticity pivot mjd',
1692 size=pars['MIRRORCHROMATICITYPIVOT'].size)
1693 parSchema.addField('compCcdChromaticity', type='ArrayD',
1694 doc='Computed CCD chromaticity terms',
1695 size=pars['COMPCCDCHROMATICITY'].size)
1696 parSchema.addField('compMedianSedSlope', type='ArrayD',
1697 doc='Computed median SED slope (per band)',
1698 size=pars['COMPMEDIANSEDSLOPE'].size)
1699 parSchema.addField('compAperCorrPivot', type='ArrayD', doc='Aperture correction pivot',
1700 size=pars['COMPAPERCORRPIVOT'].size)
1701 parSchema.addField('compAperCorrSlope', type='ArrayD', doc='Aperture correction slope',
1702 size=pars['COMPAPERCORRSLOPE'].size)
1703 parSchema.addField('compAperCorrSlopeErr', type='ArrayD', doc='Aperture correction slope error',
1704 size=pars['COMPAPERCORRSLOPEERR'].size)
1705 parSchema.addField('compAperCorrRange', type='ArrayD', doc='Aperture correction range',
1706 size=pars['COMPAPERCORRRANGE'].size)
1707 parSchema.addField('compModelErrExptimePivot', type='ArrayD', doc='Model error exptime pivot',
1708 size=pars['COMPMODELERREXPTIMEPIVOT'].size)
1709 parSchema.addField('compModelErrFwhmPivot', type='ArrayD', doc='Model error fwhm pivot',
1710 size=pars['COMPMODELERRFWHMPIVOT'].size)
1711 parSchema.addField('compModelErrSkyPivot', type='ArrayD', doc='Model error sky pivot',
1712 size=pars['COMPMODELERRSKYPIVOT'].size)
1713 parSchema.addField('compModelErrPars', type='ArrayD', doc='Model error parameters',
1714 size=pars['COMPMODELERRPARS'].size)
1715 parSchema.addField('compExpGray', type='ArrayD', doc='Computed exposure gray',
1716 size=pars['COMPEXPGRAY'].size)
1717 parSchema.addField('compVarGray', type='ArrayD', doc='Computed exposure variance',
1718 size=pars['COMPVARGRAY'].size)
1719 parSchema.addField('compExpDeltaMagBkg', type='ArrayD',
1720 doc='Computed exposure offset due to background',
1721 size=pars['COMPEXPDELTAMAGBKG'].size)
1722 parSchema.addField('compNGoodStarPerExp', type='ArrayI',
1723 doc='Computed number of good stars per exposure',
1724 size=pars['COMPNGOODSTARPEREXP'].size)
1725 parSchema.addField('compExpRefOffset', type='ArrayD',
1726 doc='Computed per-visit median offset between standard stars and ref stars.',
1727 size=pars['COMPEXPREFOFFSET'].size)
1728 parSchema.addField('compSigFgcm', type='ArrayD', doc='Computed sigma_fgcm (intrinsic repeatability)',
1729 size=pars['COMPSIGFGCM'].size)
1730 parSchema.addField('compSigmaCal', type='ArrayD', doc='Computed sigma_cal (systematic error floor)',
1731 size=pars['COMPSIGMACAL'].size)
1732 parSchema.addField('compRetrievedLnPwv', type='ArrayD', doc='Retrieved ln(pwv) (smoothed)',
1733 size=pars['COMPRETRIEVEDLNPWV'].size)
1734 parSchema.addField('compRetrievedLnPwvRaw', type='ArrayD', doc='Retrieved ln(pwv) (raw)',
1735 size=pars['COMPRETRIEVEDLNPWVRAW'].size)
1736 parSchema.addField('compRetrievedLnPwvFlag', type='ArrayI', doc='Retrieved ln(pwv) Flag',
1737 size=pars['COMPRETRIEVEDLNPWVFLAG'].size)
1738 parSchema.addField('compRetrievedTauNight', type='ArrayD', doc='Retrieved tau (per night)',
1739 size=pars['COMPRETRIEVEDTAUNIGHT'].size)
1740 parSchema.addField('compEpsilon', type='ArrayD',
1741 doc='Computed epsilon background offset per visit (nJy/arcsec2)',
1742 size=pars['COMPEPSILON'].size)
1743 parSchema.addField('compMedDeltaAper', type='ArrayD',
1744 doc='Median delta mag aper per visit',
1745 size=pars['COMPMEDDELTAAPER'].size)
1746 parSchema.addField('compGlobalEpsilon', type='ArrayD',
1747 doc='Computed epsilon bkg offset (global) (nJy/arcsec2)',
1748 size=pars['COMPGLOBALEPSILON'].size)
1749 parSchema.addField('compEpsilonMap', type='ArrayD',
1750 doc='Computed epsilon maps (nJy/arcsec2)',
1751 size=pars['COMPEPSILONMAP'].size)
1752 parSchema.addField('compEpsilonNStarMap', type='ArrayI',
1753 doc='Number of stars per pixel in computed epsilon maps',
1754 size=pars['COMPEPSILONNSTARMAP'].size)
1755 parSchema.addField('compEpsilonCcdMap', type='ArrayD',
1756 doc='Computed epsilon ccd maps (nJy/arcsec2)',
1757 size=pars['COMPEPSILONCCDMAP'].size)
1758 parSchema.addField('compEpsilonCcdNStarMap', type='ArrayI',
1759 doc='Number of stars per ccd bin in epsilon ccd maps',
1760 size=pars['COMPEPSILONCCDNSTARMAP'].size)
1761 # superstarflat section
1762 parSchema.addField('superstarSize', type='ArrayI', doc='Superstar matrix size',
1763 size=4)
1764 parSchema.addField('superstar', type='ArrayD', doc='Superstar matrix (flattened)',
1765 size=parSuperStarFlat.size)
1766
1767 return parSchema
1768
1769 def _makeParCatalog(self, parSchema, parInfo, pars, parSuperStarFlat,
1770 lutFilterNameString, fitBandString):
1771 """
1772 Make the FGCM parameter catalog for persistence
1773
1774 Parameters
1775 ----------
1776 parSchema: `lsst.afw.table.Schema`
1777 Parameter catalog schema
1778 pars: `numpy.ndarray`
1779 FGCM parameters to put into parCat
1780 parSuperStarFlat: `numpy.array`
1781 FGCM superstar flat array to put into parCat
1782 lutFilterNameString: `str`
1783 Combined string of all the lutFilterNames
1784 fitBandString: `str`
1785 Combined string of all the fitBands
1786
1787 Returns
1788 -------
1789 parCat: `afwTable.BasicCatalog`
1790 Atmosphere and instrumental model parameter catalog for persistence
1791 """
1792
1793 parCat = afwTable.BaseCatalog(parSchema)
1794 parCat.reserve(1)
1795
1796 # The parameter catalog just has one row, with many columns for all the
1797 # atmosphere and instrument fit parameters
1798 rec = parCat.addNew()
1799
1800 # info section
1801 rec['nCcd'] = parInfo['NCCD'][0]
1802 rec['lutFilterNames'] = lutFilterNameString
1803 rec['fitBands'] = fitBandString
1804 # note these are not currently supported here.
1805 rec['hasExternalPwv'] = 0
1806 rec['hasExternalTau'] = 0
1807
1808 # parameter section
1809
1810 scalarNames = ['parRetrievedLnPwvScale', 'parRetrievedLnPwvOffset']
1811
1812 arrNames = ['parAlpha', 'parO3', 'parLnTauIntercept', 'parLnTauSlope',
1813 'parLnPwvIntercept', 'parLnPwvSlope', 'parLnPwvQuadratic',
1814 'parQeSysIntercept', 'compQeSysSlope',
1815 'parRetrievedLnPwvNightlyOffset', 'compAperCorrPivot',
1816 'parFilterOffset', 'parFilterOffsetFitFlag',
1817 'compAbsThroughput', 'compRefOffset', 'compRefSigma',
1818 'compMirrorChromaticity', 'mirrorChromaticityPivot', 'compCcdChromaticity',
1819 'compAperCorrSlope', 'compAperCorrSlopeErr', 'compAperCorrRange',
1820 'compModelErrExptimePivot', 'compModelErrFwhmPivot',
1821 'compModelErrSkyPivot', 'compModelErrPars',
1822 'compExpGray', 'compVarGray', 'compNGoodStarPerExp', 'compSigFgcm',
1823 'compSigmaCal', 'compExpDeltaMagBkg', 'compMedianSedSlope',
1824 'compRetrievedLnPwv', 'compRetrievedLnPwvRaw', 'compRetrievedLnPwvFlag',
1825 'compRetrievedTauNight', 'compEpsilon', 'compMedDeltaAper',
1826 'compGlobalEpsilon', 'compEpsilonMap', 'compEpsilonNStarMap',
1827 'compEpsilonCcdMap', 'compEpsilonCcdNStarMap', 'compExpRefOffset']
1828
1829 for scalarName in scalarNames:
1830 rec[scalarName] = pars[scalarName.upper()][0]
1831
1832 for arrName in arrNames:
1833 rec[arrName][:] = np.atleast_1d(pars[0][arrName.upper()])[:]
1834
1835 # superstar section
1836 rec['superstarSize'][:] = parSuperStarFlat.shape
1837 rec['superstar'][:] = parSuperStarFlat.ravel()
1838
1839 return parCat
1840
1841 def _makeFlagStarSchema(self):
1842 """
1843 Make the flagged-stars schema
1844
1845 Returns
1846 -------
1847 flagStarSchema: `lsst.afw.table.Schema`
1848 """
1849
1850 flagStarSchema = afwTable.Schema()
1851
1852 flagStarSchema.addField('objId', type=np.int32, doc='FGCM object id')
1853 flagStarSchema.addField('objFlag', type=np.int32, doc='FGCM object flag')
1854
1855 return flagStarSchema
1856
1857 def _makeFlagStarCat(self, flagStarSchema, flagStarStruct):
1858 """
1859 Make the flagged star catalog for persistence
1860
1861 Parameters
1862 ----------
1863 flagStarSchema: `lsst.afw.table.Schema`
1864 Flagged star schema
1865 flagStarStruct: `numpy.ndarray`
1866 Flagged star structure from fgcm
1867
1868 Returns
1869 -------
1870 flagStarCat: `lsst.afw.table.BaseCatalog`
1871 Flagged star catalog for persistence
1872 """
1873
1874 flagStarCat = afwTable.BaseCatalog(flagStarSchema)
1875 flagStarCat.resize(flagStarStruct.size)
1876
1877 flagStarCat['objId'][:] = flagStarStruct['OBJID']
1878 flagStarCat['objFlag'][:] = flagStarStruct['OBJFLAG']
1879
1880 return flagStarCat
Defines the fields and offsets for a table.
Definition Schema.h:51