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