23 """Build star observations for input to FGCM.
25 This task finds all the visits and calexps in a repository (or a subset
26 based on command line parameters) and extract all the potential calibration
27 stars for input into fgcm. This task additionally uses fgcm to match
28 star observations into unique stars, and performs as much cleaning of
29 the input catalog as possible.
40 from .fgcmBuildStarsBase
import FgcmBuildStarsConfigBase, FgcmBuildStarsRunner, FgcmBuildStarsBaseTask
41 from .utilities
import computeApproxPixelAreaFields
43 __all__ = [
'FgcmBuildStarsConfig',
'FgcmBuildStarsTask']
47 """Config for FgcmBuildStarsTask"""
49 referenceCCD = pexConfig.Field(
50 doc=
"Reference CCD for scanning visits",
54 checkAllCcds = pexConfig.Field(
55 doc=(
"Check repo for all CCDs for each visit specified. To be used when the "
56 "full set of ids (visit/ccd) are not specified on the command line. For "
57 "Gen2, specifying one ccd and setting checkAllCcds=True is significantly "
58 "faster than the alternatives."),
72 fluxFlagName = self.
instFluxFieldinstFluxField[0: -len(
'instFlux')] +
'flag'
73 sourceSelector.flags.bad = [
'base_PixelFlags_flag_edge',
74 'base_PixelFlags_flag_interpolatedCenter',
75 'base_PixelFlags_flag_saturatedCenter',
76 'base_PixelFlags_flag_crCenter',
77 'base_PixelFlags_flag_bad',
78 'base_PixelFlags_flag_interpolated',
79 'base_PixelFlags_flag_saturated',
85 sourceSelector.flags.bad.append(localBackgroundFlagName)
87 sourceSelector.signalToNoise.fluxField = self.
instFluxFieldinstFluxField
88 sourceSelector.signalToNoise.errField = self.
instFluxFieldinstFluxField +
'Err'
93 Build stars for the FGCM global calibration, using src catalogs.
95 ConfigClass = FgcmBuildStarsConfig
96 RunnerClass = FgcmBuildStarsRunner
97 _DefaultName =
"fgcmBuildStars"
100 def _makeArgumentParser(cls):
101 """Create an argument parser"""
102 parser = pipeBase.ArgumentParser(name=cls.
_DefaultName_DefaultName)
103 parser.add_id_argument(
"--id",
"src", help=
"Data ID, e.g. --id visit=6789")
107 def _findAndGroupDataRefsGen2(self, butler, camera, dataRefs):
108 self.log.
info(
"Grouping dataRefs by %s" % (self.config.visitDataRefName))
111 for detector
in camera:
112 ccdIds.append(detector.getId())
121 for dataRef
in dataRefs:
122 visit = dataRef.dataId[self.config.visitDataRefName]
124 if not dataRef.datasetExists(datasetType=
'src'):
127 if self.config.checkAllCcds:
128 if visit
in groupedDataRefs:
131 dataId = dataRef.dataId.copy()
134 dataId[self.config.ccdDataRefName] = ccdId
135 if butler.datasetExists(
'src', dataId=dataId):
136 goodDataRef = butler.dataRef(
'src', dataId=dataId)
137 if visit
in groupedDataRefs:
138 if (goodDataRef.dataId[self.config.ccdDataRefName]
not in
139 [d.dataId[self.config.ccdDataRefName]
for d
in groupedDataRefs[visit]]):
140 groupedDataRefs[visit].
append(goodDataRef)
144 groupedDataRefs[visit] = [goodDataRef]
148 if visit
in groupedDataRefs:
149 if (dataRef.dataId[self.config.ccdDataRefName]
not in
150 [d.dataId[self.config.ccdDataRefName]
for d
in groupedDataRefs[visit]]):
151 groupedDataRefs[visit].
append(dataRef)
155 groupedDataRefs[visit] = [dataRef]
157 if (nVisits % 100) == 0
and nVisits > 0:
158 self.log.
info(
"Found %d unique %ss..." % (nVisits,
159 self.config.visitDataRefName))
161 self.log.
info(
"Found %d unique %ss total." % (nVisits,
162 self.config.visitDataRefName))
165 def ccdSorter(dataRef):
166 ccdId = dataRef.dataId[self.config.ccdDataRefName]
167 if ccdId == self.config.referenceCCD:
173 if not self.config.checkAllCcds:
174 for visit
in groupedDataRefs:
175 groupedDataRefs[visit] = sorted(groupedDataRefs[visit], key=ccdSorter)
178 return dict(sorted(groupedDataRefs.items()))
183 calibFluxApertureRadius=None,
184 visitCatDataRef=None,
187 startTime = time.time()
193 if (visitCatDataRef
is not None and starObsDataRef
is None
194 or visitCatDataRef
is None and starObsDataRef
is not None):
195 self.log.
warning(
"Only one of visitCatDataRef and starObsDataRef are set, so "
196 "no checkpoint files will be persisted.")
198 if self.config.doSubtractLocalBackground
and calibFluxApertureRadius
is None:
199 raise RuntimeError(
"Must set calibFluxApertureRadius if doSubtractLocalBackground is True.")
202 dataRef = groupedDataRefs[
list(groupedDataRefs.keys())[0]][0]
206 for ccdIndex, detector
in enumerate(camera):
207 ccdMapping[detector.getId()] = ccdIndex
216 outputSchema = sourceMapper.getOutputSchema()
218 if inStarObsCat
is not None:
219 fullCatalog = inStarObsCat
220 comp1 = fullCatalog.schema.compare(outputSchema, outputSchema.EQUAL_KEYS)
221 comp2 = fullCatalog.schema.compare(outputSchema, outputSchema.EQUAL_NAMES)
222 if not comp1
or not comp2:
223 raise RuntimeError(
"Existing fgcmStarObservations file found with mismatched schema.")
229 instFluxKey = sourceSchema[self.config.instFluxField].asKey()
230 instFluxErrKey = sourceSchema[self.config.instFluxField +
'Err'].asKey()
231 visitKey = outputSchema[
'visit'].asKey()
232 ccdKey = outputSchema[
'ccd'].asKey()
233 instMagKey = outputSchema[
'instMag'].asKey()
234 instMagErrKey = outputSchema[
'instMagErr'].asKey()
235 deltaMagBkgKey = outputSchema[
'deltaMagBkg'].asKey()
238 if self.config.doSubtractLocalBackground:
239 localBackgroundFluxKey = sourceSchema[self.config.localBackgroundFluxField].asKey()
240 localBackgroundArea = np.pi*calibFluxApertureRadius**2.
242 aperOutputSchema = aperMapper.getOutputSchema()
244 instFluxAperInKey = sourceSchema[self.config.apertureInnerInstFluxField].asKey()
245 instFluxErrAperInKey = sourceSchema[self.config.apertureInnerInstFluxField +
'Err'].asKey()
246 instFluxAperOutKey = sourceSchema[self.config.apertureOuterInstFluxField].asKey()
247 instFluxErrAperOutKey = sourceSchema[self.config.apertureOuterInstFluxField +
'Err'].asKey()
248 instMagInKey = aperOutputSchema[
'instMag_aper_inner'].asKey()
249 instMagErrInKey = aperOutputSchema[
'instMagErr_aper_inner'].asKey()
250 instMagOutKey = aperOutputSchema[
'instMag_aper_outer'].asKey()
251 instMagErrOutKey = aperOutputSchema[
'instMagErr_aper_outer'].asKey()
256 for ctr, visit
in enumerate(visitCat):
257 if visit[
'sources_read']:
260 expTime = visit[
'exptime']
267 for dataRef
in groupedDataRefs[visit[
'visit']]:
269 ccdId = dataRef.dataId[self.config.ccdDataRefName]
271 sources = dataRef.get(datasetType=
'src', flags=afwTable.SOURCE_IO_NO_FOOTPRINTS)
272 goodSrc = self.sourceSelector.selectSources(sources)
276 if self.config.doSubtractLocalBackground:
277 localBackground = localBackgroundArea*sources[localBackgroundFluxKey]
279 bad, = np.where((sources[instFluxKey] - localBackground) <= 0.0)
280 goodSrc.selected[bad] =
False
283 tempCat.reserve(goodSrc.selected.sum())
284 tempCat.extend(sources[goodSrc.selected], mapper=sourceMapper)
285 tempCat[visitKey][:] = visit[
'visit']
286 tempCat[ccdKey][:] = ccdId
289 scaledInstFlux = (sources[instFluxKey][goodSrc.selected]
290 * visit[
'scaling'][ccdMapping[ccdId]])
291 tempCat[instMagKey][:] = (-2.5*np.log10(scaledInstFlux) + 2.5*np.log10(expTime))
294 if self.config.doSubtractLocalBackground:
308 tempCat[deltaMagBkgKey][:] = (-2.5*np.log10(sources[instFluxKey][goodSrc.selected]
309 - localBackground[goodSrc.selected]) -
310 -2.5*np.log10(sources[instFluxKey][goodSrc.selected]))
312 tempCat[deltaMagBkgKey][:] = 0.0
317 tempCat[instMagErrKey][:] = k*(sources[instFluxErrKey][goodSrc.selected]
318 / sources[instFluxKey][goodSrc.selected])
321 tempCat[
'jacobian'] = approxPixelAreaFields[ccdId].evaluate(tempCat[
'x'],
325 if self.config.doApplyWcsJacobian:
326 tempCat[instMagKey][:] -= 2.5*np.log10(tempCat[
'jacobian'][:])
328 fullCatalog.extend(tempCat)
333 tempAperCat.reserve(goodSrc.selected.sum())
334 tempAperCat.extend(sources[goodSrc.selected], mapper=aperMapper)
336 with np.warnings.catch_warnings():
339 np.warnings.simplefilter(
"ignore")
341 tempAperCat[instMagInKey][:] = -2.5*np.log10(
342 sources[instFluxAperInKey][goodSrc.selected])
343 tempAperCat[instMagErrInKey][:] = k*(
344 sources[instFluxErrAperInKey][goodSrc.selected]
345 / sources[instFluxAperInKey][goodSrc.selected])
346 tempAperCat[instMagOutKey][:] = -2.5*np.log10(
347 sources[instFluxAperOutKey][goodSrc.selected])
348 tempAperCat[instMagErrOutKey][:] = k*(
349 sources[instFluxErrAperOutKey][goodSrc.selected]
350 / sources[instFluxAperOutKey][goodSrc.selected])
352 aperVisitCatalog.extend(tempAperCat)
354 nStarInVisit += len(tempCat)
357 if not aperVisitCatalog.isContiguous():
358 aperVisitCatalog = aperVisitCatalog.copy(deep=
True)
360 instMagIn = aperVisitCatalog[instMagInKey]
361 instMagErrIn = aperVisitCatalog[instMagErrInKey]
362 instMagOut = aperVisitCatalog[instMagOutKey]
363 instMagErrOut = aperVisitCatalog[instMagErrOutKey]
365 ok = (np.isfinite(instMagIn) & np.isfinite(instMagErrIn)
366 & np.isfinite(instMagOut) & np.isfinite(instMagErrOut))
368 visit[
'deltaAper'] = np.median(instMagIn[ok] - instMagOut[ok])
369 visit[
'sources_read'] =
True
371 self.log.
info(
" Found %d good stars in visit %d (deltaAper = %.3f)" %
372 (nStarInVisit, visit[
'visit'], visit[
'deltaAper']))
374 if ((ctr % self.config.nVisitsPerCheckpoint) == 0
375 and starObsDataRef
is not None and visitCatDataRef
is not None):
378 starObsDataRef.put(fullCatalog)
379 visitCatDataRef.put(visitCat)
381 self.log.
info(
"Found all good star observations in %.2f s" %
382 (time.time() - startTime))
386 def _makeAperMapper(self, sourceSchema):
388 Make a schema mapper for fgcm aperture measurements
392 sourceSchema: `afwTable.Schema`
393 Default source schema from the butler
397 aperMapper: `afwTable.schemaMapper`
398 Mapper to the FGCM aperture schema
402 aperMapper.addMapping(sourceSchema[
'coord_ra'].asKey(),
'ra')
403 aperMapper.addMapping(sourceSchema[
'coord_dec'].asKey(),
'dec')
404 aperMapper.editOutputSchema().addField(
'instMag_aper_inner', type=np.float64,
405 doc=
"Magnitude at inner aperture")
406 aperMapper.editOutputSchema().addField(
'instMagErr_aper_inner', type=np.float64,
407 doc=
"Magnitude error at inner aperture")
408 aperMapper.editOutputSchema().addField(
'instMag_aper_outer', type=np.float64,
409 doc=
"Magnitude at outer aperture")
410 aperMapper.editOutputSchema().addField(
'instMagErr_aper_outer', type=np.float64,
411 doc=
"Magnitude error at outer aperture")
A mapping between the keys of two Schemas, used to copy data between them.
def _makeAperMapper(self, sourceSchema)
def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat, sourceSchema, camera, calibFluxApertureRadius=None, visitCatDataRef=None, starObsDataRef=None, inStarObsCat=None)
def _makeSourceMapper(self, sourceSchema)
doSubtractLocalBackground
daf::base::PropertyList * list
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
def computeApproxPixelAreaFields(camera)