LSST Applications g013ef56533+7c9321ec0f,g042eb84c57+c6cfa41bc3,g199a45376c+0ba108daf9,g1fd858c14a+fcad0d0313,g210f2d0738+c0f94c6586,g262e1987ae+a7e710680e,g29ae962dfc+fb55f2edb0,g2ac17093b6+61d6563b1e,g2b1d02342f+df6f932764,g2cef7863aa+aef1011c0b,g2f7ad74990+c0f94c6586,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53cf87ae69,g47891489e3+4316d04fff,g511e8cfd20+baa56acf6c,g53246c7159+8c5ae1fdc5,g54cd7ddccb+fd7ad03fde,g64539dfbff+c0f94c6586,g67b6fd64d1+4316d04fff,g67fd3c3899+c0f94c6586,g6985122a63+4316d04fff,g74acd417e5+ca833bee28,g786e29fd12+668abc6043,g81db2e9a8d+b2ec8e584f,g87389fa792+8856018cbb,g89139ef638+4316d04fff,g8d7436a09f+0a24083b20,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+11eb8fd5b8,gbf99507273+8c5ae1fdc5,gcdda8b9158+e4c84c9d5c,gce8aa8abaa+8c5ae1fdc5,gd7ef33dd92+4316d04fff,gdab6d2f7ff+ca833bee28,ge410e46f29+4316d04fff,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.40
LSST Data Management Base Package
Loading...
Searching...
No Matches
fgcmCalibrateTractBase.py
Go to the documentation of this file.
1# This file is part of fgcmcal.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21"""Base class for running fgcmcal on a single tract using src tables
22or sourceTable_visit tables.
23"""
24
25import abc
26
27import numpy as np
28
29import lsst.pex.config as pexConfig
30import lsst.pipe.base as pipeBase
31
32from .fgcmBuildStarsTable import FgcmBuildStarsTableTask
33from .fgcmFitCycle import FgcmFitCycleConfig
34from .fgcmOutputProducts import FgcmOutputProductsTask
35from .utilities import makeConfigDict, translateFgcmLut, translateVisitCatalog
36from .utilities import computeApertureRadiusFromName, extractReferenceMags
37from .utilities import makeZptSchema, makeZptCat
38from .utilities import makeAtmSchema, makeAtmCat
39from .utilities import makeStdSchema, makeStdCat
40from .focalPlaneProjector import FocalPlaneProjector
41
42import fgcm
43
44__all__ = ['FgcmCalibrateTractConfigBase', 'FgcmCalibrateTractBaseTask']
45
46
47class FgcmCalibrateTractConfigBase(pexConfig.Config):
48 """Config for FgcmCalibrateTract"""
49
50 fgcmBuildStars = pexConfig.ConfigurableField(
51 target=FgcmBuildStarsTableTask,
52 doc="Task to load and match stars for fgcm",
53 )
54 fgcmFitCycle = pexConfig.ConfigField(
55 dtype=FgcmFitCycleConfig,
56 doc="Config to run a single fgcm fit cycle",
57 )
58 fgcmOutputProducts = pexConfig.ConfigurableField(
59 target=FgcmOutputProductsTask,
60 doc="Task to output fgcm products",
61 )
62 convergenceTolerance = pexConfig.Field(
63 doc="Tolerance on repeatability convergence (per band)",
64 dtype=float,
65 default=0.005,
66 )
67 maxFitCycles = pexConfig.Field(
68 doc="Maximum number of fit cycles",
69 dtype=int,
70 default=5,
71 )
72 doDebuggingPlots = pexConfig.Field(
73 doc="Make plots for debugging purposes?",
74 dtype=bool,
75 default=False,
76 )
77
78 def setDefaults(self):
79 pexConfig.Config.setDefaults(self)
80
81 self.fgcmFitCycle.quietMode = True
82 self.fgcmFitCycle.doPlots = False
83 self.fgcmOutputProducts.doReferenceCalibration = False
84 self.fgcmOutputProducts.photoCal.applyColorTerms = False
85
86 def validate(self):
87 super().validate()
88
89 for band in self.fgcmFitCycle.bands:
90 if not self.fgcmFitCycle.useRepeatabilityForExpGrayCutsDict[band]:
91 msg = 'Must set useRepeatabilityForExpGrayCutsDict[band]=True for all bands'
92 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict,
93 self, msg)
94
95
96class FgcmCalibrateTractBaseTask(pipeBase.PipelineTask, abc.ABC):
97 """Base class to calibrate a single tract using fgcmcal
98 """
99 def __init__(self, initInputs=None, **kwargs):
100 super().__init__(**kwargs)
101 self.makeSubtask("fgcmBuildStars", initInputs=initInputs)
102 self.makeSubtask("fgcmOutputProducts")
103
104 def run(self, handleDict, tract,
105 buildStarsRefObjLoader=None, returnCatalogs=True):
106 """Run the calibrations for a single tract with fgcm.
107
108 Parameters
109 ----------
110 handleDict : `dict`
111 All handles are `lsst.daf.butler.DeferredDatasetHandle`
112 handle dictionary with the following keys. Note that all
113 keys need not be set based on config parameters.
114
115 ``"camera"``
116 Camera object (`lsst.afw.cameraGeom.Camera`)
117 ``"source_catalogs"``
118 `list` of handles for input source catalogs.
119 ``"sourceSchema"``
120 Schema for the source catalogs.
121 ``"fgcmLookUpTable"``
122 handle for the FGCM look-up table.
123 ``"calexps"``
124 `list` of handles for the input calexps
125 ``"fgcmPhotoCalibs"``
126 `dict` of output photoCalib handles. Key is
127 (tract, visit, detector).
128 Present if doZeropointOutput is True.
129 ``"fgcmTransmissionAtmospheres"``
130 `dict` of output atmosphere transmission handles.
131 Key is (tract, visit).
132 Present if doAtmosphereOutput is True.
133 tract : `int`
134 Tract number
135 buildStarsRefObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`, optional
136 Reference object loader object for fgcmBuildStars.
137 returnCatalogs : `bool`, optional
138 Return photoCalibs as per-visit exposure catalogs.
139
140 Returns
141 -------
142 outstruct : `lsst.pipe.base.Struct`
143 Output structure with keys:
144
145 offsets : `np.ndarray`
146 Final reference offsets, per band.
147 repeatability : `np.ndarray`
148 Raw fgcm repeatability for bright stars, per band.
149 atmospheres : `generator` [(`int`, `lsst.afw.image.TransmissionCurve`)]
150 Generator that returns (visit, transmissionCurve) tuples.
151 photoCalibs : `generator` [(`int`, `int`, `str`, `lsst.afw.image.PhotoCalib`)]
152 Generator that returns (visit, ccd, filtername, photoCalib) tuples.
153 (returned if returnCatalogs is False).
154 photoCalibCatalogs : `generator` [(`int`, `lsst.afw.table.ExposureCatalog`)]
155 Generator that returns (visit, exposureCatalog) tuples.
156 (returned if returnCatalogs is True).
157 """
158 self.log.info("Running on tract %d", (tract))
159
160 # Compute the aperture radius if necessary. This is useful to do now before
161 # any heavy lifting has happened (fail early).
162 calibFluxApertureRadius = None
163 if self.config.fgcmBuildStars.doSubtractLocalBackground:
164 try:
165 field = self.config.fgcmBuildStars.instFluxField
166 calibFluxApertureRadius = computeApertureRadiusFromName(field)
167 except RuntimeError:
168 raise RuntimeError("Could not determine aperture radius from %s. "
169 "Cannot use doSubtractLocalBackground." %
170 (field))
171
172 # Run the build stars tasks
173
174 # Note that we will need visitCat at the end of the procedure for the outputs
175 groupedHandles = self.fgcmBuildStars._groupHandles(handleDict['sourceTableHandleDict'],
176 handleDict['visitSummaryHandleDict'])
177
178 lutCat = handleDict["fgcmLookUpTable"].get()
179 camera = handleDict["camera"]
180 if len(camera) == lutCat[0]["nCcd"]:
181 useScienceDetectors = False
182 else:
183 # If the LUT has a different number of detectors than
184 # the camera, then we only want to use science detectors
185 # in the focal plane projector.
186 useScienceDetectors = True
187 del lutCat
188
189 visitCat = self.fgcmBuildStars.fgcmMakeVisitCatalog(
190 camera,
191 groupedHandles,
192 useScienceDetectors=useScienceDetectors,
193 )
194 rad = calibFluxApertureRadius
195 fgcmStarObservationCat = self.fgcmBuildStars.fgcmMakeAllStarObservations(groupedHandles,
196 visitCat,
197 handleDict['sourceSchema'],
198 handleDict['camera'],
199 calibFluxApertureRadius=rad)
200
201 if self.fgcmBuildStars.config.doReferenceMatches:
202 lutHandle = handleDict['fgcmLookUpTable']
203 self.fgcmBuildStars.makeSubtask("fgcmLoadReferenceCatalog",
204 refObjLoader=buildStarsRefObjLoader,
205 refCatName=self.fgcmBuildStars.config.connections.refCat)
206 else:
207 lutHandle = None
208
209 fgcmStarIdCat, fgcmStarIndicesCat, fgcmRefCat = \
210 self.fgcmBuildStars.fgcmMatchStars(visitCat,
211 fgcmStarObservationCat,
212 lutHandle=lutHandle)
213
214 # Load the LUT
215 lutCat = handleDict['fgcmLookUpTable'].get()
216 fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat,
217 dict(self.config.fgcmFitCycle.physicalFilterMap))
218 del lutCat
219
220 # Translate the visit catalog into fgcm format
221 fgcmExpInfo = translateVisitCatalog(visitCat)
222
223 configDict = makeConfigDict(self.config.fgcmFitCycle, self.log, handleDict['camera'],
224 self.config.fgcmFitCycle.maxIterBeforeFinalCycle,
225 True, False, lutIndexVals[0]['FILTERNAMES'],
226 tract=tract)
227
228 focalPlaneProjector = FocalPlaneProjector(handleDict['camera'],
229 self.config.fgcmFitCycle.defaultCameraOrientation)
230
231 # Set up the fit cycle task
232
233 noFitsDict = {'lutIndex': lutIndexVals,
234 'lutStd': lutStd,
235 'expInfo': fgcmExpInfo,
236 'focalPlaneProjector': focalPlaneProjector}
237
238 fgcmFitCycle = fgcm.FgcmFitCycle(configDict, useFits=False,
239 noFitsDict=noFitsDict, noOutput=True)
240
241 # We determine the conversion from the native units (typically radians) to
242 # degrees for the first star. This allows us to treat coord_ra/coord_dec as
243 # numpy arrays rather than Angles, which would we approximately 600x slower.
244 conv = fgcmStarObservationCat[0]['ra'].asDegrees() / float(fgcmStarObservationCat[0]['ra'])
245
246 # To load the stars, we need an initial parameter object
247 fgcmPars = fgcm.FgcmParameters.newParsWithArrays(fgcmFitCycle.fgcmConfig,
248 fgcmLut,
249 fgcmExpInfo)
250
251 # Match star observations to visits
252 # Only those star observations that match visits from fgcmExpInfo['VISIT'] will
253 # actually be transferred into fgcm using the indexing below.
254
255 obsIndex = fgcmStarIndicesCat['obsIndex']
256 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'],
257 fgcmStarObservationCat['visit'][obsIndex])
258
259 refMag, refMagErr = extractReferenceMags(fgcmRefCat,
260 self.config.fgcmFitCycle.bands,
261 self.config.fgcmFitCycle.physicalFilterMap)
262 refId = fgcmRefCat['fgcm_id'][:]
263
264 fgcmStars = fgcm.FgcmStars(fgcmFitCycle.fgcmConfig)
265 fgcmStars.loadStars(fgcmPars,
266 fgcmStarObservationCat['visit'][obsIndex],
267 fgcmStarObservationCat['ccd'][obsIndex],
268 fgcmStarObservationCat['ra'][obsIndex] * conv,
269 fgcmStarObservationCat['dec'][obsIndex] * conv,
270 fgcmStarObservationCat['instMag'][obsIndex],
271 fgcmStarObservationCat['instMagErr'][obsIndex],
272 fgcmExpInfo['FILTERNAME'][visitIndex],
273 fgcmStarIdCat['fgcm_id'][:],
274 fgcmStarIdCat['ra'][:],
275 fgcmStarIdCat['dec'][:],
276 fgcmStarIdCat['obsArrIndex'][:],
277 fgcmStarIdCat['nObs'][:],
278 obsX=fgcmStarObservationCat['x'][obsIndex],
279 obsY=fgcmStarObservationCat['y'][obsIndex],
280 obsDeltaMagBkg=fgcmStarObservationCat['deltaMagBkg'][obsIndex],
281 obsDeltaAper=fgcmStarObservationCat['deltaMagAper'][obsIndex],
282 psfCandidate=fgcmStarObservationCat['psf_candidate'][obsIndex],
283 refID=refId,
284 refMag=refMag,
285 refMagErr=refMagErr,
286 flagID=None,
287 flagFlag=None,
288 computeNobs=True)
289
290 # Clear out some memory
291 del fgcmStarIdCat
292 del fgcmStarIndicesCat
293 del fgcmRefCat
294
295 fgcmFitCycle.setLUT(fgcmLut)
296 fgcmFitCycle.setStars(fgcmStars, fgcmPars)
297
298 converged = False
299 cycleNumber = 0
300
301 previousReservedRawRepeatability = np.zeros(fgcmPars.nBands) + 1000.0
302 previousParInfo = None
303 previousParams = None
304 previousSuperStar = None
305
306 while (not converged and cycleNumber < self.config.maxFitCycles):
307
308 fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
309
310 if cycleNumber > 0:
311 # Use parameters from previous cycle
312 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
313 fgcmExpInfo,
314 previousParInfo,
315 previousParams,
316 previousSuperStar)
317 # We need to reset the star magnitudes and errors for the next
318 # cycle
319 fgcmFitCycle.fgcmStars.reloadStarMagnitudes(fgcmStarObservationCat['instMag'][obsIndex],
320 fgcmStarObservationCat['instMagErr'][obsIndex])
321 fgcmFitCycle.initialCycle = False
322
323 fgcmFitCycle.setPars(fgcmPars)
324 fgcmFitCycle.finishSetup()
325
326 fgcmFitCycle.run()
327
328 # Grab the parameters for the next cycle
329 previousParInfo, previousParams = fgcmFitCycle.fgcmPars.parsToArrays()
330 previousSuperStar = fgcmFitCycle.fgcmPars.parSuperStarFlat.copy()
331
332 self.log.info("Raw repeatability after cycle number %d is:" % (cycleNumber))
333 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
334 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
335 continue
336 rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
337 self.log.info(" Band %s, repeatability: %.2f mmag" % (band, rep))
338
339 # Check for convergence
340 if np.all((previousReservedRawRepeatability
341 - fgcmFitCycle.fgcmPars.compReservedRawRepeatability)
342 < self.config.convergenceTolerance):
343 self.log.info("Raw repeatability has converged after cycle number %d." % (cycleNumber))
344 converged = True
345 else:
346 fgcmFitCycle.fgcmConfig.expGrayPhotometricCut[:] = fgcmFitCycle.updatedPhotometricCut
347 fgcmFitCycle.fgcmConfig.expGrayHighCut[:] = fgcmFitCycle.updatedHighCut
348 fgcmFitCycle.fgcmConfig.precomputeSuperStarInitialCycle = False
349 fgcmFitCycle.fgcmConfig.freezeStdAtmosphere = False
350 previousReservedRawRepeatability[:] = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
351 self.log.info("Setting exposure gray photometricity cuts to:")
352 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
353 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
354 continue
355 cut = fgcmFitCycle.updatedPhotometricCut[i] * 1000.0
356 self.log.info(" Band %s, photometricity cut: %.2f mmag" % (band, cut))
357
358 cycleNumber += 1
359
360 # Log warning if not converged
361 if not converged:
362 self.log.warning("Maximum number of fit cycles exceeded (%d) without convergence.", cycleNumber)
363
364 # Do final clean-up iteration
365 fgcmFitCycle.fgcmConfig.freezeStdAtmosphere = False
366 fgcmFitCycle.fgcmConfig.resetParameters = False
367 fgcmFitCycle.fgcmConfig.maxIter = 0
368 fgcmFitCycle.fgcmConfig.outputZeropoints = True
369 fgcmFitCycle.fgcmConfig.outputStandards = True
370 fgcmFitCycle.fgcmConfig.doPlots = self.config.doDebuggingPlots
371 fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
372 fgcmFitCycle.initialCycle = False
373
374 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
375 fgcmExpInfo,
376 previousParInfo,
377 previousParams,
378 previousSuperStar)
379 fgcmFitCycle.fgcmStars.reloadStarMagnitudes(fgcmStarObservationCat['instMag'][obsIndex],
380 fgcmStarObservationCat['instMagErr'][obsIndex])
381 fgcmFitCycle.setPars(fgcmPars)
382 fgcmFitCycle.finishSetup()
383
384 self.log.info("Running final clean-up fit cycle...")
385 fgcmFitCycle.run()
386
387 self.log.info("Raw repeatability after clean-up cycle is:")
388 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
389 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
390 continue
391 rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
392 self.log.info(" Band %s, repeatability: %.2f mmag" % (band, rep))
393
394 # Do the outputs. Need to keep track of tract.
395
396 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
397 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
398
399 zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
400 zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
401
402 atmSchema = makeAtmSchema()
403 atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
404
405 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
406 stdSchema = makeStdSchema(len(goodBands))
407 stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
408
409 outStruct = self.fgcmOutputProducts.generateTractOutputProducts(handleDict,
410 tract,
411 visitCat,
412 zptCat, atmCat, stdCat,
413 self.config.fgcmBuildStars)
414
415 outStruct.repeatability = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
416
417 fgcmFitCycle.freeSharedMemory()
418
419 return outStruct
run(self, handleDict, tract, buildStarsRefObjLoader=None, returnCatalogs=True)