LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
fgcmBuildStars.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 """Build star observations for input to FGCM.
24 
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.
30 """
31 
32 import time
33 
34 import numpy as np
35 
36 import lsst.pex.config as pexConfig
37 import lsst.pipe.base as pipeBase
38 import lsst.afw.table as afwTable
39 
40 from .fgcmBuildStarsBase import FgcmBuildStarsConfigBase, FgcmBuildStarsRunner, FgcmBuildStarsBaseTask
41 from .utilities import computeApproxPixelAreaFields
42 
43 __all__ = ['FgcmBuildStarsConfig', 'FgcmBuildStarsTask']
44 
45 
47  """Config for FgcmBuildStarsTask"""
48 
49  referenceCCD = pexConfig.Field(
50  doc="Reference CCD for scanning visits",
51  dtype=int,
52  default=13,
53  )
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."),
59  dtype=bool,
60  default=True,
61  )
62 
63  def setDefaults(self):
64  super().setDefaults()
65 
66  sourceSelector = self.sourceSelectorsourceSelector["science"]
67 
68  # The names here correspond to raw src catalogs, which differ
69  # from the post-transformed sourceTable_visit catalogs.
70  # Therefore, field and flag names cannot be easily
71  # derived from the base config class.
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',
80  'slot_Centroid_flag',
81  fluxFlagName]
82 
83  if self.doSubtractLocalBackgrounddoSubtractLocalBackground:
84  localBackgroundFlagName = self.localBackgroundFluxFieldlocalBackgroundFluxField[0: -len('instFlux')] + 'flag'
85  sourceSelector.flags.bad.append(localBackgroundFlagName)
86 
87  sourceSelector.signalToNoise.fluxField = self.instFluxFieldinstFluxField
88  sourceSelector.signalToNoise.errField = self.instFluxFieldinstFluxField + 'Err'
89 
90 
92  """
93  Build stars for the FGCM global calibration, using src catalogs.
94  """
95  ConfigClass = FgcmBuildStarsConfig
96  RunnerClass = FgcmBuildStarsRunner
97  _DefaultName = "fgcmBuildStars"
98 
99  @classmethod
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")
104 
105  return parser
106 
107  def _findAndGroupDataRefsGen2(self, butler, camera, dataRefs):
108  self.log.info("Grouping dataRefs by %s" % (self.config.visitDataRefName))
109 
110  ccdIds = []
111  for detector in camera:
112  ccdIds.append(detector.getId())
113 
114  # TODO: related to DM-13730, this dance of looking for source visits
115  # will be unnecessary with Gen3 Butler. This should be part of
116  # DM-13730.
117 
118  nVisits = 0
119 
120  groupedDataRefs = {}
121  for dataRef in dataRefs:
122  visit = dataRef.dataId[self.config.visitDataRefName]
123  # If we don't have the dataset, just continue
124  if not dataRef.datasetExists(datasetType='src'):
125  continue
126  # If we need to check all ccds, do it here
127  if self.config.checkAllCcds:
128  if visit in groupedDataRefs:
129  # We already have found this visit
130  continue
131  dataId = dataRef.dataId.copy()
132  # For each ccd we must check that a valid source catalog exists.
133  for ccdId in ccdIds:
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)
141  else:
142  # This is a new visit
143  nVisits += 1
144  groupedDataRefs[visit] = [goodDataRef]
145  else:
146  # We have already confirmed that the dataset exists, so no need
147  # to check here.
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)
152  else:
153  # This is a new visit
154  nVisits += 1
155  groupedDataRefs[visit] = [dataRef]
156 
157  if (nVisits % 100) == 0 and nVisits > 0:
158  self.log.info("Found %d unique %ss..." % (nVisits,
159  self.config.visitDataRefName))
160 
161  self.log.info("Found %d unique %ss total." % (nVisits,
162  self.config.visitDataRefName))
163 
164  # Put them in ccd order, with the reference ccd first (if available)
165  def ccdSorter(dataRef):
166  ccdId = dataRef.dataId[self.config.ccdDataRefName]
167  if ccdId == self.config.referenceCCD:
168  return -100
169  else:
170  return ccdId
171 
172  # If we did not check all ccds, put them in ccd order
173  if not self.config.checkAllCcds:
174  for visit in groupedDataRefs:
175  groupedDataRefs[visit] = sorted(groupedDataRefs[visit], key=ccdSorter)
176 
177  # This should be sorted by visit (the key)
178  return dict(sorted(groupedDataRefs.items()))
179 
180  def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat,
181  sourceSchema,
182  camera,
183  calibFluxApertureRadius=None,
184  visitCatDataRef=None,
185  starObsDataRef=None,
186  inStarObsCat=None):
187  startTime = time.time()
188 
189  # If both dataRefs are None, then we assume the caller does not
190  # want to store checkpoint files. If both are set, we will
191  # do checkpoint files. And if only one is set, this is potentially
192  # unintentional and we will warn.
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.")
197 
198  if self.config.doSubtractLocalBackground and calibFluxApertureRadius is None:
199  raise RuntimeError("Must set calibFluxApertureRadius if doSubtractLocalBackground is True.")
200 
201  # create our source schema. Use the first valid dataRef
202  dataRef = groupedDataRefs[list(groupedDataRefs.keys())[0]][0]
203 
204  # Construct a mapping from ccd number to index
205  ccdMapping = {}
206  for ccdIndex, detector in enumerate(camera):
207  ccdMapping[detector.getId()] = ccdIndex
208 
209  approxPixelAreaFields = computeApproxPixelAreaFields(camera)
210 
211  sourceMapper = self._makeSourceMapper_makeSourceMapper(sourceSchema)
212 
213  # We also have a temporary catalog that will accumulate aperture measurements
214  aperMapper = self._makeAperMapper_makeAperMapper(sourceSchema)
215 
216  outputSchema = sourceMapper.getOutputSchema()
217 
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.")
224  else:
225  fullCatalog = afwTable.BaseCatalog(outputSchema)
226 
227  # FGCM will provide relative calibration for the flux in config.instFluxField
228 
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()
236 
237  # Prepare local background if desired
238  if self.config.doSubtractLocalBackground:
239  localBackgroundFluxKey = sourceSchema[self.config.localBackgroundFluxField].asKey()
240  localBackgroundArea = np.pi*calibFluxApertureRadius**2.
241 
242  aperOutputSchema = aperMapper.getOutputSchema()
243 
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()
252 
253  k = 2.5/np.log(10.)
254 
255  # loop over visits
256  for ctr, visit in enumerate(visitCat):
257  if visit['sources_read']:
258  continue
259 
260  expTime = visit['exptime']
261 
262  nStarInVisit = 0
263 
264  # Reset the aperture catalog (per visit)
265  aperVisitCatalog = afwTable.BaseCatalog(aperOutputSchema)
266 
267  for dataRef in groupedDataRefs[visit['visit']]:
268 
269  ccdId = dataRef.dataId[self.config.ccdDataRefName]
270 
271  sources = dataRef.get(datasetType='src', flags=afwTable.SOURCE_IO_NO_FOOTPRINTS)
272  goodSrc = self.sourceSelector.selectSources(sources)
273 
274  # Need to add a selection based on the local background correction
275  # if necessary
276  if self.config.doSubtractLocalBackground:
277  localBackground = localBackgroundArea*sources[localBackgroundFluxKey]
278 
279  bad, = np.where((sources[instFluxKey] - localBackground) <= 0.0)
280  goodSrc.selected[bad] = False
281 
282  tempCat = afwTable.BaseCatalog(fullCatalog.schema)
283  tempCat.reserve(goodSrc.selected.sum())
284  tempCat.extend(sources[goodSrc.selected], mapper=sourceMapper)
285  tempCat[visitKey][:] = visit['visit']
286  tempCat[ccdKey][:] = ccdId
287 
288  # Compute "instrumental magnitude" by scaling flux with exposure time.
289  scaledInstFlux = (sources[instFluxKey][goodSrc.selected]
290  * visit['scaling'][ccdMapping[ccdId]])
291  tempCat[instMagKey][:] = (-2.5*np.log10(scaledInstFlux) + 2.5*np.log10(expTime))
292 
293  # Compute the change in magnitude from the background offset
294  if self.config.doSubtractLocalBackground:
295  # At the moment we only adjust the flux and not the flux
296  # error by the background because the error on
297  # base_LocalBackground_instFlux is the rms error in the
298  # background annulus, not the error on the mean in the
299  # background estimate (which is much smaller, by sqrt(n)
300  # pixels used to estimate the background, which we do not
301  # have access to in this task). In the default settings,
302  # the annulus is sufficiently large such that these
303  # additional errors are are negligibly small (much less
304  # than a mmag in quadrature).
305 
306  # This is the difference between the mag with background correction
307  # and the mag without background correction.
308  tempCat[deltaMagBkgKey][:] = (-2.5*np.log10(sources[instFluxKey][goodSrc.selected]
309  - localBackground[goodSrc.selected]) -
310  -2.5*np.log10(sources[instFluxKey][goodSrc.selected]))
311  else:
312  tempCat[deltaMagBkgKey][:] = 0.0
313 
314  # Compute instMagErr from instFluxErr/instFlux, any scaling
315  # will cancel out.
316 
317  tempCat[instMagErrKey][:] = k*(sources[instFluxErrKey][goodSrc.selected]
318  / sources[instFluxKey][goodSrc.selected])
319 
320  # Compute the jacobian from an approximate PixelAreaBoundedField
321  tempCat['jacobian'] = approxPixelAreaFields[ccdId].evaluate(tempCat['x'],
322  tempCat['y'])
323 
324  # Apply the jacobian if configured
325  if self.config.doApplyWcsJacobian:
326  tempCat[instMagKey][:] -= 2.5*np.log10(tempCat['jacobian'][:])
327 
328  fullCatalog.extend(tempCat)
329 
330  # And the aperture information
331  # This does not need the jacobian because it is all locally relative
332  tempAperCat = afwTable.BaseCatalog(aperVisitCatalog.schema)
333  tempAperCat.reserve(goodSrc.selected.sum())
334  tempAperCat.extend(sources[goodSrc.selected], mapper=aperMapper)
335 
336  with np.warnings.catch_warnings():
337  # Ignore warnings, we will filter infinities and
338  # nans below.
339  np.warnings.simplefilter("ignore")
340 
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])
351 
352  aperVisitCatalog.extend(tempAperCat)
353 
354  nStarInVisit += len(tempCat)
355 
356  # Compute the median delta-aper
357  if not aperVisitCatalog.isContiguous():
358  aperVisitCatalog = aperVisitCatalog.copy(deep=True)
359 
360  instMagIn = aperVisitCatalog[instMagInKey]
361  instMagErrIn = aperVisitCatalog[instMagErrInKey]
362  instMagOut = aperVisitCatalog[instMagOutKey]
363  instMagErrOut = aperVisitCatalog[instMagErrOutKey]
364 
365  ok = (np.isfinite(instMagIn) & np.isfinite(instMagErrIn)
366  & np.isfinite(instMagOut) & np.isfinite(instMagErrOut))
367 
368  visit['deltaAper'] = np.median(instMagIn[ok] - instMagOut[ok])
369  visit['sources_read'] = True
370 
371  self.log.info(" Found %d good stars in visit %d (deltaAper = %.3f)" %
372  (nStarInVisit, visit['visit'], visit['deltaAper']))
373 
374  if ((ctr % self.config.nVisitsPerCheckpoint) == 0
375  and starObsDataRef is not None and visitCatDataRef is not None):
376  # We need to persist both the stars and the visit catalog which gets
377  # additional metadata from each visit.
378  starObsDataRef.put(fullCatalog)
379  visitCatDataRef.put(visitCat)
380 
381  self.log.info("Found all good star observations in %.2f s" %
382  (time.time() - startTime))
383 
384  return fullCatalog
385 
386  def _makeAperMapper(self, sourceSchema):
387  """
388  Make a schema mapper for fgcm aperture measurements
389 
390  Parameters
391  ----------
392  sourceSchema: `afwTable.Schema`
393  Default source schema from the butler
394 
395  Returns
396  -------
397  aperMapper: `afwTable.schemaMapper`
398  Mapper to the FGCM aperture schema
399  """
400 
401  aperMapper = afwTable.SchemaMapper(sourceSchema)
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")
412 
413  return aperMapper
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat, sourceSchema, camera, calibFluxApertureRadius=None, visitCatDataRef=None, starObsDataRef=None, inStarObsCat=None)
daf::base::PropertyList * list
Definition: fits.cc:913
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
def computeApproxPixelAreaFields(camera)
Definition: utilities.py:499