LSST Applications g1635faa6d4+215bc75b8c,g1653933729+a8ce1bb630,g22ce9dc20b+d972d8df89,g28da252d5a+0fcf840c6d,g29321ee8c0+e558be0e74,g2bbee38e9b+9634bc57db,g2bc492864f+9634bc57db,g2cdde0e794+c2c89b37c4,g3156d2b45e+41e33cbcdc,g347aa1857d+9634bc57db,g35bb328faa+a8ce1bb630,g3a166c0a6a+9634bc57db,g3e281a1b8c+9f2c4e2fc3,g414038480c+077ccc18e7,g41af890bb2+2a6f257a1d,g5fbc88fb19+17cd334064,g781aacb6e4+a8ce1bb630,g7ab3e175f3+59ce30aec6,g80478fca09+f8b2ab54e1,g82479be7b0+ba9d578ff8,g858d7b2824+d972d8df89,g9125e01d80+a8ce1bb630,g9726552aa6+10f999ec6a,ga5288a1d22+630363936d,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gb9c6c11c1e+9553554aa7,gbd46683f8f+0c4209622a,gc28159a63d+9634bc57db,gcf0d15dbbd+2db122af0a,gda3e153d99+d972d8df89,gda6a2b7d83+2db122af0a,gdaeeff99f8+1711a396fd,ge2409df99d+d1dc2f3b25,ge33fd446bb+d972d8df89,ge79ae78c31+9634bc57db,gf0baf85859+147a0692ba,gf3967379c6+02b11634a5,w.2024.45
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 visitCat = self.fgcmBuildStars.fgcmMakeVisitCatalog(handleDict['camera'], groupedHandles)
178 rad = calibFluxApertureRadius
179 fgcmStarObservationCat = self.fgcmBuildStars.fgcmMakeAllStarObservations(groupedHandles,
180 visitCat,
181 handleDict['sourceSchema'],
182 handleDict['camera'],
183 calibFluxApertureRadius=rad)
184
185 if self.fgcmBuildStars.config.doReferenceMatches:
186 lutHandle = handleDict['fgcmLookUpTable']
187 self.fgcmBuildStars.makeSubtask("fgcmLoadReferenceCatalog",
188 refObjLoader=buildStarsRefObjLoader,
189 refCatName=self.fgcmBuildStars.config.connections.refCat)
190 else:
191 lutHandle = None
192
193 fgcmStarIdCat, fgcmStarIndicesCat, fgcmRefCat = \
194 self.fgcmBuildStars.fgcmMatchStars(visitCat,
195 fgcmStarObservationCat,
196 lutHandle=lutHandle)
197
198 # Load the LUT
199 lutCat = handleDict['fgcmLookUpTable'].get()
200 fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat,
201 dict(self.config.fgcmFitCycle.physicalFilterMap))
202 del lutCat
203
204 # Translate the visit catalog into fgcm format
205 fgcmExpInfo = translateVisitCatalog(visitCat)
206
207 configDict = makeConfigDict(self.config.fgcmFitCycle, self.log, handleDict['camera'],
208 self.config.fgcmFitCycle.maxIterBeforeFinalCycle,
209 True, False, lutIndexVals[0]['FILTERNAMES'],
210 tract=tract)
211
212 focalPlaneProjector = FocalPlaneProjector(handleDict['camera'],
213 self.config.fgcmFitCycle.defaultCameraOrientation)
214
215 # Set up the fit cycle task
216
217 noFitsDict = {'lutIndex': lutIndexVals,
218 'lutStd': lutStd,
219 'expInfo': fgcmExpInfo,
220 'focalPlaneProjector': focalPlaneProjector}
221
222 fgcmFitCycle = fgcm.FgcmFitCycle(configDict, useFits=False,
223 noFitsDict=noFitsDict, noOutput=True)
224
225 # We determine the conversion from the native units (typically radians) to
226 # degrees for the first star. This allows us to treat coord_ra/coord_dec as
227 # numpy arrays rather than Angles, which would we approximately 600x slower.
228 conv = fgcmStarObservationCat[0]['ra'].asDegrees() / float(fgcmStarObservationCat[0]['ra'])
229
230 # To load the stars, we need an initial parameter object
231 fgcmPars = fgcm.FgcmParameters.newParsWithArrays(fgcmFitCycle.fgcmConfig,
232 fgcmLut,
233 fgcmExpInfo)
234
235 # Match star observations to visits
236 # Only those star observations that match visits from fgcmExpInfo['VISIT'] will
237 # actually be transferred into fgcm using the indexing below.
238
239 obsIndex = fgcmStarIndicesCat['obsIndex']
240 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'],
241 fgcmStarObservationCat['visit'][obsIndex])
242
243 refMag, refMagErr = extractReferenceMags(fgcmRefCat,
244 self.config.fgcmFitCycle.bands,
245 self.config.fgcmFitCycle.physicalFilterMap)
246 refId = fgcmRefCat['fgcm_id'][:]
247
248 fgcmStars = fgcm.FgcmStars(fgcmFitCycle.fgcmConfig)
249 fgcmStars.loadStars(fgcmPars,
250 fgcmStarObservationCat['visit'][obsIndex],
251 fgcmStarObservationCat['ccd'][obsIndex],
252 fgcmStarObservationCat['ra'][obsIndex] * conv,
253 fgcmStarObservationCat['dec'][obsIndex] * conv,
254 fgcmStarObservationCat['instMag'][obsIndex],
255 fgcmStarObservationCat['instMagErr'][obsIndex],
256 fgcmExpInfo['FILTERNAME'][visitIndex],
257 fgcmStarIdCat['fgcm_id'][:],
258 fgcmStarIdCat['ra'][:],
259 fgcmStarIdCat['dec'][:],
260 fgcmStarIdCat['obsArrIndex'][:],
261 fgcmStarIdCat['nObs'][:],
262 obsX=fgcmStarObservationCat['x'][obsIndex],
263 obsY=fgcmStarObservationCat['y'][obsIndex],
264 obsDeltaMagBkg=fgcmStarObservationCat['deltaMagBkg'][obsIndex],
265 obsDeltaAper=fgcmStarObservationCat['deltaMagAper'][obsIndex],
266 psfCandidate=fgcmStarObservationCat['psf_candidate'][obsIndex],
267 refID=refId,
268 refMag=refMag,
269 refMagErr=refMagErr,
270 flagID=None,
271 flagFlag=None,
272 computeNobs=True)
273
274 # Clear out some memory
275 del fgcmStarIdCat
276 del fgcmStarIndicesCat
277 del fgcmRefCat
278
279 fgcmFitCycle.setLUT(fgcmLut)
280 fgcmFitCycle.setStars(fgcmStars, fgcmPars)
281
282 converged = False
283 cycleNumber = 0
284
285 previousReservedRawRepeatability = np.zeros(fgcmPars.nBands) + 1000.0
286 previousParInfo = None
287 previousParams = None
288 previousSuperStar = None
289
290 while (not converged and cycleNumber < self.config.maxFitCycles):
291
292 fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
293
294 if cycleNumber > 0:
295 # Use parameters from previous cycle
296 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
297 fgcmExpInfo,
298 previousParInfo,
299 previousParams,
300 previousSuperStar)
301 # We need to reset the star magnitudes and errors for the next
302 # cycle
303 fgcmFitCycle.fgcmStars.reloadStarMagnitudes(fgcmStarObservationCat['instMag'][obsIndex],
304 fgcmStarObservationCat['instMagErr'][obsIndex])
305 fgcmFitCycle.initialCycle = False
306
307 fgcmFitCycle.setPars(fgcmPars)
308 fgcmFitCycle.finishSetup()
309
310 fgcmFitCycle.run()
311
312 # Grab the parameters for the next cycle
313 previousParInfo, previousParams = fgcmFitCycle.fgcmPars.parsToArrays()
314 previousSuperStar = fgcmFitCycle.fgcmPars.parSuperStarFlat.copy()
315
316 self.log.info("Raw repeatability after cycle number %d is:" % (cycleNumber))
317 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
318 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
319 continue
320 rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
321 self.log.info(" Band %s, repeatability: %.2f mmag" % (band, rep))
322
323 # Check for convergence
324 if np.all((previousReservedRawRepeatability
325 - fgcmFitCycle.fgcmPars.compReservedRawRepeatability)
326 < self.config.convergenceTolerance):
327 self.log.info("Raw repeatability has converged after cycle number %d." % (cycleNumber))
328 converged = True
329 else:
330 fgcmFitCycle.fgcmConfig.expGrayPhotometricCut[:] = fgcmFitCycle.updatedPhotometricCut
331 fgcmFitCycle.fgcmConfig.expGrayHighCut[:] = fgcmFitCycle.updatedHighCut
332 fgcmFitCycle.fgcmConfig.precomputeSuperStarInitialCycle = False
333 fgcmFitCycle.fgcmConfig.freezeStdAtmosphere = False
334 previousReservedRawRepeatability[:] = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
335 self.log.info("Setting exposure gray photometricity cuts to:")
336 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
337 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
338 continue
339 cut = fgcmFitCycle.updatedPhotometricCut[i] * 1000.0
340 self.log.info(" Band %s, photometricity cut: %.2f mmag" % (band, cut))
341
342 cycleNumber += 1
343
344 # Log warning if not converged
345 if not converged:
346 self.log.warning("Maximum number of fit cycles exceeded (%d) without convergence.", cycleNumber)
347
348 # Do final clean-up iteration
349 fgcmFitCycle.fgcmConfig.freezeStdAtmosphere = False
350 fgcmFitCycle.fgcmConfig.resetParameters = False
351 fgcmFitCycle.fgcmConfig.maxIter = 0
352 fgcmFitCycle.fgcmConfig.outputZeropoints = True
353 fgcmFitCycle.fgcmConfig.outputStandards = True
354 fgcmFitCycle.fgcmConfig.doPlots = self.config.doDebuggingPlots
355 fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
356 fgcmFitCycle.initialCycle = False
357
358 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
359 fgcmExpInfo,
360 previousParInfo,
361 previousParams,
362 previousSuperStar)
363 fgcmFitCycle.fgcmStars.reloadStarMagnitudes(fgcmStarObservationCat['instMag'][obsIndex],
364 fgcmStarObservationCat['instMagErr'][obsIndex])
365 fgcmFitCycle.setPars(fgcmPars)
366 fgcmFitCycle.finishSetup()
367
368 self.log.info("Running final clean-up fit cycle...")
369 fgcmFitCycle.run()
370
371 self.log.info("Raw repeatability after clean-up cycle is:")
372 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
373 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
374 continue
375 rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
376 self.log.info(" Band %s, repeatability: %.2f mmag" % (band, rep))
377
378 # Do the outputs. Need to keep track of tract.
379
380 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
381 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
382
383 zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
384 zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
385
386 atmSchema = makeAtmSchema()
387 atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
388
389 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
390 stdSchema = makeStdSchema(len(goodBands))
391 stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
392
393 outStruct = self.fgcmOutputProducts.generateTractOutputProducts(handleDict,
394 tract,
395 visitCat,
396 zptCat, atmCat, stdCat,
397 self.config.fgcmBuildStars)
398
399 outStruct.repeatability = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
400
401 fgcmFitCycle.freeSharedMemory()
402
403 return outStruct
run(self, handleDict, tract, buildStarsRefObjLoader=None, returnCatalogs=True)