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
fgcmBuildStarsBase.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 BuildStars using src tables or sourceTable_visit tables.
22 """
23 
24 import os
25 import sys
26 import traceback
27 import abc
28 
29 import numpy as np
30 
31 import lsst.daf.persistence as dafPersist
32 import lsst.pex.config as pexConfig
33 import lsst.pipe.base as pipeBase
34 import lsst.afw.table as afwTable
35 import lsst.geom as geom
36 from lsst.daf.base import PropertyList
37 from lsst.daf.base.dateTime import DateTime
38 from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry
39 from lsst.utils.timer import timeMethod
40 
41 from .utilities import computeApertureRadiusFromDataRef
42 from .fgcmLoadReferenceCatalog import FgcmLoadReferenceCatalogTask
43 
44 import fgcm
45 
46 REFSTARS_FORMAT_VERSION = 1
47 
48 __all__ = ['FgcmBuildStarsConfigBase', 'FgcmBuildStarsRunner', 'FgcmBuildStarsBaseTask']
49 
50 
51 class FgcmBuildStarsConfigBase(pexConfig.Config):
52  """Base config for FgcmBuildStars tasks"""
53 
54  instFluxField = pexConfig.Field(
55  doc=("Faull name of the source instFlux field to use, including 'instFlux'. "
56  "The associated flag will be implicitly included in badFlags"),
57  dtype=str,
58  default='slot_CalibFlux_instFlux',
59  )
60  minPerBand = pexConfig.Field(
61  doc="Minimum observations per band",
62  dtype=int,
63  default=2,
64  )
65  matchRadius = pexConfig.Field(
66  doc="Match radius (arcseconds)",
67  dtype=float,
68  default=1.0,
69  )
70  isolationRadius = pexConfig.Field(
71  doc="Isolation radius (arcseconds)",
72  dtype=float,
73  default=2.0,
74  )
75  densityCutNside = pexConfig.Field(
76  doc="Density cut healpix nside",
77  dtype=int,
78  default=128,
79  )
80  densityCutMaxPerPixel = pexConfig.Field(
81  doc="Density cut number of stars per pixel",
82  dtype=int,
83  default=1000,
84  )
85  randomSeed = pexConfig.Field(
86  doc="Random seed for high density down-sampling.",
87  dtype=int,
88  default=None,
89  optional=True,
90  )
91  matchNside = pexConfig.Field(
92  doc="Healpix Nside for matching",
93  dtype=int,
94  default=4096,
95  )
96  coarseNside = pexConfig.Field(
97  doc="Healpix coarse Nside for partitioning matches",
98  dtype=int,
99  default=8,
100  )
101  # The following config will not be necessary after Gen2 retirement.
102  # In the meantime, obs packages should set to 'filterDefinitions.filter_to_band'
103  # which is easiest to access in the config file.
104  physicalFilterMap = pexConfig.DictField(
105  doc="Mapping from 'physicalFilter' to band.",
106  keytype=str,
107  itemtype=str,
108  default={},
109  )
110  requiredBands = pexConfig.ListField(
111  doc="Bands required for each star",
112  dtype=str,
113  default=(),
114  )
115  primaryBands = pexConfig.ListField(
116  doc=("Bands for 'primary' star matches. "
117  "A star must be observed in one of these bands to be considered "
118  "as a calibration star."),
119  dtype=str,
120  default=None
121  )
122  visitDataRefName = pexConfig.Field(
123  doc="dataRef name for the 'visit' field, usually 'visit'.",
124  dtype=str,
125  default="visit"
126  )
127  ccdDataRefName = pexConfig.Field(
128  doc="dataRef name for the 'ccd' field, usually 'ccd' or 'detector'.",
129  dtype=str,
130  default="ccd"
131  )
132  doApplyWcsJacobian = pexConfig.Field(
133  doc="Apply the jacobian of the WCS to the star observations prior to fit?",
134  dtype=bool,
135  default=True
136  )
137  doModelErrorsWithBackground = pexConfig.Field(
138  doc="Model flux errors with background term?",
139  dtype=bool,
140  default=True
141  )
142  psfCandidateName = pexConfig.Field(
143  doc="Name of field with psf candidate flag for propagation",
144  dtype=str,
145  default="calib_psf_candidate"
146  )
147  doSubtractLocalBackground = pexConfig.Field(
148  doc=("Subtract the local background before performing calibration? "
149  "This is only supported for circular aperture calibration fluxes."),
150  dtype=bool,
151  default=False
152  )
153  localBackgroundFluxField = pexConfig.Field(
154  doc="Full name of the local background instFlux field to use.",
155  dtype=str,
156  default='base_LocalBackground_instFlux'
157  )
158  sourceSelector = sourceSelectorRegistry.makeField(
159  doc="How to select sources",
160  default="science"
161  )
162  apertureInnerInstFluxField = pexConfig.Field(
163  doc=("Full name of instFlux field that contains inner aperture "
164  "flux for aperture correction proxy"),
165  dtype=str,
166  default='base_CircularApertureFlux_12_0_instFlux'
167  )
168  apertureOuterInstFluxField = pexConfig.Field(
169  doc=("Full name of instFlux field that contains outer aperture "
170  "flux for aperture correction proxy"),
171  dtype=str,
172  default='base_CircularApertureFlux_17_0_instFlux'
173  )
174  doReferenceMatches = pexConfig.Field(
175  doc="Match reference catalog as additional constraint on calibration",
176  dtype=bool,
177  default=True,
178  )
179  fgcmLoadReferenceCatalog = pexConfig.ConfigurableField(
180  target=FgcmLoadReferenceCatalogTask,
181  doc="FGCM reference object loader",
182  )
183  nVisitsPerCheckpoint = pexConfig.Field(
184  doc="Number of visits read between checkpoints",
185  dtype=int,
186  default=500,
187  )
188 
189  def setDefaults(self):
190  sourceSelector = self.sourceSelectorsourceSelector["science"]
191  sourceSelector.setDefaults()
192 
193  sourceSelector.doFlags = True
194  sourceSelector.doUnresolved = True
195  sourceSelector.doSignalToNoise = True
196  sourceSelector.doIsolated = True
197 
198  sourceSelector.signalToNoise.minimum = 10.0
199  sourceSelector.signalToNoise.maximum = 1000.0
200 
201  # FGCM operates on unresolved sources, and this setting is
202  # appropriate for the current base_ClassificationExtendedness
203  sourceSelector.unresolved.maximum = 0.5
204 
205 
206 class FgcmBuildStarsRunner(pipeBase.ButlerInitializedTaskRunner):
207  """Subclass of TaskRunner for FgcmBuildStars tasks
208 
209  fgcmBuildStarsTask.run() and fgcmBuildStarsTableTask.run() take a number of
210  arguments, one of which is the butler (for persistence and mapper data),
211  and a list of dataRefs extracted from the command line. Note that FGCM
212  runs on a large set of dataRefs, and not on single dataRef/tract/patch.
213  This class transforms the process arguments generated by the ArgumentParser
214  into the arguments expected by FgcmBuildStarsTask.run(). This runner does
215  not use any parallelization.
216  """
217  @staticmethod
218  def getTargetList(parsedCmd):
219  """
220  Return a list with one element: a tuple with the butler and
221  list of dataRefs
222  """
223  # we want to combine the butler with any (or no!) dataRefs
224  return [(parsedCmd.butler, parsedCmd.id.refList)]
225 
226  def __call__(self, args):
227  """
228  Parameters
229  ----------
230  args: `tuple` with (butler, dataRefList)
231 
232  Returns
233  -------
234  exitStatus: `list` with `lsst.pipe.base.Struct`
235  exitStatus (0: success; 1: failure)
236  """
237  butler, dataRefList = args
238 
239  task = self.TaskClass(config=self.config, log=self.log)
240 
241  exitStatus = 0
242  if self.doRaise:
243  task.runDataRef(butler, dataRefList)
244  else:
245  try:
246  task.runDataRef(butler, dataRefList)
247  except Exception as e:
248  exitStatus = 1
249  task.log.fatal("Failed: %s" % e)
250  if not isinstance(e, pipeBase.TaskError):
251  traceback.print_exc(file=sys.stderr)
252 
253  task.writeMetadata(butler)
254 
255  # The task does not return any results:
256  return [pipeBase.Struct(exitStatus=exitStatus)]
257 
258  def run(self, parsedCmd):
259  """
260  Run the task, with no multiprocessing
261 
262  Parameters
263  ----------
264  parsedCmd: `lsst.pipe.base.ArgumentParser` parsed command line
265  """
266 
267  resultList = []
268 
269  if self.precall(parsedCmd):
270  targetList = self.getTargetListgetTargetList(parsedCmd)
271  resultList = self(targetList[0])
272 
273  return resultList
274 
275 
276 class FgcmBuildStarsBaseTask(pipeBase.PipelineTask, pipeBase.CmdLineTask, abc.ABC):
277  """
278  Base task to build stars for FGCM global calibration
279 
280  Parameters
281  ----------
282  butler : `lsst.daf.persistence.Butler`
283  """
284  def __init__(self, initInputs=None, butler=None, **kwargs):
285  super().__init__(**kwargs)
286 
287  self.makeSubtask("sourceSelector")
288  # Only log warning and fatal errors from the sourceSelector
289  self.sourceSelector.log.setLevel(self.sourceSelector.log.WARN)
290 
291  # no saving of metadata for now
292  def _getMetadataName(self):
293  return None
294 
295  @timeMethod
296  def runDataRef(self, butler, dataRefs):
297  """
298  Cross-match and make star list for FGCM Input
299 
300  Parameters
301  ----------
302  butler: `lsst.daf.persistence.Butler`
303  dataRefs: `list` of `lsst.daf.persistence.ButlerDataRef`
304  Source data references for the input visits.
305 
306  Raises
307  ------
308  RuntimeErrror: Raised if `config.doReferenceMatches` is set and
309  an fgcmLookUpTable is not available, or if computeFluxApertureRadius()
310  fails if the calibFlux is not a CircularAperture flux.
311  """
312  datasetType = dataRefs[0].butlerSubset.datasetType
313  self.log.info("Running with %d %s dataRefs", len(dataRefs), datasetType)
314 
315  if self.config.doReferenceMatches:
316  self.makeSubtask("fgcmLoadReferenceCatalog", butler=butler)
317  # Ensure that we have a LUT
318  if not butler.datasetExists('fgcmLookUpTable'):
319  raise RuntimeError("Must have fgcmLookUpTable if using config.doReferenceMatches")
320  # Compute aperture radius if necessary. This is useful to do now before
321  # any heavy lifting has happened (fail early).
322  calibFluxApertureRadius = None
323  if self.config.doSubtractLocalBackground:
324  try:
325  calibFluxApertureRadius = computeApertureRadiusFromDataRef(dataRefs[0],
326  self.config.instFluxField)
327  except RuntimeError as e:
328  raise RuntimeError("Could not determine aperture radius from %s. "
329  "Cannot use doSubtractLocalBackground." %
330  (self.config.instFluxField)) from e
331 
332  camera = butler.get('camera')
333  groupedDataRefs = self._findAndGroupDataRefsGen2_findAndGroupDataRefsGen2(butler, camera, dataRefs)
334 
335  # Make the visit catalog if necessary
336  # First check if the visit catalog is in the _current_ path
337  # We cannot use Gen2 datasetExists() because that checks all parent
338  # directories as well, which would make recovering from faults
339  # and fgcmcal reruns impossible.
340  visitCatDataRef = butler.dataRef('fgcmVisitCatalog')
341  filename = visitCatDataRef.get('fgcmVisitCatalog_filename')[0]
342  if os.path.exists(filename):
343  # This file exists and we should continue processing
344  inVisitCat = visitCatDataRef.get()
345  if len(inVisitCat) != len(groupedDataRefs):
346  raise RuntimeError("Existing visitCatalog found, but has an inconsistent "
347  "number of visits. Cannot continue.")
348  else:
349  inVisitCat = None
350 
351  visitCat = self.fgcmMakeVisitCatalogfgcmMakeVisitCatalog(camera, groupedDataRefs,
352  visitCatDataRef=visitCatDataRef,
353  inVisitCat=inVisitCat)
354 
355  # Persist the visitCat as a checkpoint file.
356  visitCatDataRef.put(visitCat)
357 
358  starObsDataRef = butler.dataRef('fgcmStarObservations')
359  filename = starObsDataRef.get('fgcmStarObservations_filename')[0]
360  if os.path.exists(filename):
361  inStarObsCat = starObsDataRef.get()
362  else:
363  inStarObsCat = None
364 
365  rad = calibFluxApertureRadius
366  sourceSchemaDataRef = butler.dataRef('src_schema')
367  sourceSchema = sourceSchemaDataRef.get('src_schema', immediate=True).schema
368  fgcmStarObservationCat = self.fgcmMakeAllStarObservationsfgcmMakeAllStarObservations(groupedDataRefs,
369  visitCat,
370  sourceSchema,
371  camera,
372  calibFluxApertureRadius=rad,
373  starObsDataRef=starObsDataRef,
374  visitCatDataRef=visitCatDataRef,
375  inStarObsCat=inStarObsCat)
376  visitCatDataRef.put(visitCat)
377  starObsDataRef.put(fgcmStarObservationCat)
378 
379  # Always do the matching.
380  if self.config.doReferenceMatches:
381  lutDataRef = butler.dataRef('fgcmLookUpTable')
382  else:
383  lutDataRef = None
384  fgcmStarIdCat, fgcmStarIndicesCat, fgcmRefCat = self.fgcmMatchStarsfgcmMatchStars(visitCat,
385  fgcmStarObservationCat,
386  lutDataRef=lutDataRef)
387 
388  # Persist catalogs via the butler
389  butler.put(fgcmStarIdCat, 'fgcmStarIds')
390  butler.put(fgcmStarIndicesCat, 'fgcmStarIndices')
391  if fgcmRefCat is not None:
392  butler.put(fgcmRefCat, 'fgcmReferenceStars')
393 
394  @abc.abstractmethod
395  def _findAndGroupDataRefsGen2(self, butler, camera, dataRefs):
396  """
397  Find and group dataRefs (by visit); Gen2 only.
398 
399  Parameters
400  ----------
401  butler : `lsst.daf.persistence.Butler`
402  Gen2 butler.
403  camera : `lsst.afw.cameraGeom.Camera`
404  Camera from the butler.
405  dataRefs : `list` of `lsst.daf.persistence.ButlerDataRef`
406  Data references for the input visits.
407 
408  Returns
409  -------
410  groupedDataRefs : `dict` [`int`, `list`]
411  Dictionary with sorted visit keys, and `list`s of
412  `lsst.daf.persistence.ButlerDataRef`
413  """
414  raise NotImplementedError("_findAndGroupDataRefsGen2 not implemented.")
415 
416  @abc.abstractmethod
417  def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat,
418  sourceSchema,
419  camera,
420  calibFluxApertureRadius=None,
421  visitCatDataRef=None,
422  starObsDataRef=None,
423  inStarObsCat=None):
424  """
425  Compile all good star observations from visits in visitCat. Checkpoint files
426  will be stored if both visitCatDataRef and starObsDataRef are not None.
427 
428  Parameters
429  ----------
430  groupedDataRefs : `dict` of `list`s
431  Lists of `~lsst.daf.persistence.ButlerDataRef` or
432  `~lsst.daf.butler.DeferredDatasetHandle`, grouped by visit.
433  visitCat : `~afw.table.BaseCatalog`
434  Catalog with visit data for FGCM
435  sourceSchema : `~lsst.afw.table.Schema`
436  Schema for the input src catalogs.
437  camera : `~lsst.afw.cameraGeom.Camera`
438  calibFluxApertureRadius : `float`, optional
439  Aperture radius for calibration flux.
440  visitCatDataRef : `~lsst.daf.persistence.ButlerDataRef`, optional
441  Dataref to write visitCat for checkpoints
442  starObsDataRef : `~lsst.daf.persistence.ButlerDataRef`, optional
443  Dataref to write the star observation catalog for checkpoints.
444  inStarObsCat : `~afw.table.BaseCatalog`
445  Input observation catalog. If this is incomplete, observations
446  will be appended from when it was cut off.
447 
448  Returns
449  -------
450  fgcmStarObservations : `afw.table.BaseCatalog`
451  Full catalog of good observations.
452 
453  Raises
454  ------
455  RuntimeError: Raised if doSubtractLocalBackground is True and
456  calibFluxApertureRadius is not set.
457  """
458  raise NotImplementedError("fgcmMakeAllStarObservations not implemented.")
459 
460  def fgcmMakeVisitCatalog(self, camera, groupedDataRefs, bkgDataRefDict=None,
461  visitCatDataRef=None, inVisitCat=None):
462  """
463  Make a visit catalog with all the keys from each visit
464 
465  Parameters
466  ----------
467  camera: `lsst.afw.cameraGeom.Camera`
468  Camera from the butler
469  groupedDataRefs: `dict`
470  Dictionary with visit keys, and `list`s of
471  `lsst.daf.persistence.ButlerDataRef`
472  bkgDataRefDict: `dict`, optional
473  Dictionary of gen3 dataRefHandles for background info.
474  visitCatDataRef: `lsst.daf.persistence.ButlerDataRef`, optional
475  Dataref to write visitCat for checkpoints
476  inVisitCat: `afw.table.BaseCatalog`, optional
477  Input (possibly incomplete) visit catalog
478 
479  Returns
480  -------
481  visitCat: `afw.table.BaseCatalog`
482  """
483 
484  self.log.info("Assembling visitCatalog from %d %ss" %
485  (len(groupedDataRefs), self.config.visitDataRefName))
486 
487  nCcd = len(camera)
488 
489  if inVisitCat is None:
490  schema = self._makeFgcmVisitSchema_makeFgcmVisitSchema(nCcd)
491 
492  visitCat = afwTable.BaseCatalog(schema)
493  visitCat.reserve(len(groupedDataRefs))
494  visitCat.resize(len(groupedDataRefs))
495 
496  visitCat['visit'] = list(groupedDataRefs.keys())
497  visitCat['used'] = 0
498  visitCat['sources_read'] = False
499  else:
500  visitCat = inVisitCat
501 
502  # No matter what, fill the catalog. This will check if it was
503  # already read.
504  self._fillVisitCatalog_fillVisitCatalog(visitCat, groupedDataRefs,
505  bkgDataRefDict=bkgDataRefDict,
506  visitCatDataRef=visitCatDataRef)
507 
508  return visitCat
509 
510  def _fillVisitCatalog(self, visitCat, groupedDataRefs, bkgDataRefDict=None,
511  visitCatDataRef=None):
512  """
513  Fill the visit catalog with visit metadata
514 
515  Parameters
516  ----------
517  visitCat : `afw.table.BaseCatalog`
518  Visit catalog. See _makeFgcmVisitSchema() for schema definition.
519  groupedDataRefs : `dict`
520  Dictionary with visit keys, and `list`s of
521  `lsst.daf.persistence.ButlerDataRef` or
522  `lsst.daf.butler.DeferredDatasetHandle`
523  visitCatDataRef : `lsst.daf.persistence.ButlerDataRef`, optional
524  Dataref to write ``visitCat`` for checkpoints. Gen2 only.
525  bkgDataRefDict : `dict`, optional
526  Dictionary of Gen3 `lsst.daf.butler.DeferredDatasetHandle`
527  for background info.
528  """
529  bbox = geom.BoxI(geom.PointI(0, 0), geom.PointI(1, 1))
530 
531  for i, visit in enumerate(groupedDataRefs):
532  # We don't use the bypasses since we need the psf info which does
533  # not have a bypass
534  # TODO: When DM-15500 is implemented in the Gen3 Butler, this
535  # can be fixed
536 
537  # Do not read those that have already been read
538  if visitCat['used'][i]:
539  continue
540 
541  if (i % self.config.nVisitsPerCheckpoint) == 0:
542  self.log.info("Retrieving metadata for %s %d (%d/%d)" %
543  (self.config.visitDataRefName, visit, i, len(groupedDataRefs)))
544  # Save checkpoint if desired
545  if visitCatDataRef is not None:
546  visitCatDataRef.put(visitCat)
547 
548  dataRef = groupedDataRefs[visit][0]
549  if isinstance(dataRef, dafPersist.ButlerDataRef):
550  # Gen2: calexp dataRef
551  # The first dataRef in the group will be the reference ccd (if available)
552  exp = dataRef.get(datasetType='calexp_sub', bbox=bbox)
553  visitInfo = exp.getInfo().getVisitInfo()
554  label = dataRef.get(datasetType='calexp_filterLabel')
555  physicalFilter = label.physicalLabel
556  psf = exp.getPsf()
557  bbox = dataRef.get(datasetType='calexp_bbox')
558  psfSigma = psf.computeShape(bbox.getCenter()).getDeterminantRadius()
559  else:
560  # Gen3: use the visitSummary dataRef
561  summary = dataRef.get()
562 
563  summaryRow = summary.find(self.config.referenceCCD)
564  if summaryRow is None:
565  # Take the first available ccd if reference isn't available
566  summaryRow = summary[0]
567 
568  summaryDetector = summaryRow['id']
569  visitInfo = summaryRow.getVisitInfo()
570  physicalFilter = summaryRow['physical_filter']
571  # Compute the median psf sigma if possible
572  goodSigma, = np.where(summary['psfSigma'] > 0)
573  if goodSigma.size > 2:
574  psfSigma = np.median(summary['psfSigma'][goodSigma])
575  elif goodSigma > 0:
576  psfSigma = np.mean(summary['psfSigma'][goodSigma])
577  else:
578  psfSigma = 0.0
579 
580  rec = visitCat[i]
581  rec['visit'] = visit
582  rec['physicalFilter'] = physicalFilter
583  # TODO DM-26991: when gen2 is removed, gen3 workflow will make it
584  # much easier to get the wcs's necessary to recompute the pointing
585  # ra/dec at the center of the camera.
586  radec = visitInfo.getBoresightRaDec()
587  rec['telra'] = radec.getRa().asDegrees()
588  rec['teldec'] = radec.getDec().asDegrees()
589  rec['telha'] = visitInfo.getBoresightHourAngle().asDegrees()
590  rec['telrot'] = visitInfo.getBoresightRotAngle().asDegrees()
591  rec['mjd'] = visitInfo.getDate().get(system=DateTime.MJD)
592  rec['exptime'] = visitInfo.getExposureTime()
593  # convert from Pa to millibar
594  # Note that I don't know if this unit will need to be per-camera config
595  rec['pmb'] = visitInfo.getWeather().getAirPressure() / 100
596  # Flag to signify if this is a "deep" field. Not currently used
597  rec['deepFlag'] = 0
598  # Relative flat scaling (1.0 means no relative scaling)
599  rec['scaling'][:] = 1.0
600  # Median delta aperture, to be measured from stars
601  rec['deltaAper'] = 0.0
602  rec['psfSigma'] = psfSigma
603 
604  if self.config.doModelErrorsWithBackground:
605  foundBkg = False
606  if isinstance(dataRef, dafPersist.ButlerDataRef):
607  # Gen2-style dataRef
608  det = dataRef.dataId[self.config.ccdDataRefName]
609  if dataRef.datasetExists(datasetType='calexpBackground'):
610  bgList = dataRef.get(datasetType='calexpBackground')
611  foundBkg = True
612  else:
613  # Gen3-style dataRef
614  try:
615  # Use the same detector used from the summary.
616  bkgRef = bkgDataRefDict[(visit, summaryDetector)]
617  bgList = bkgRef.get()
618  foundBkg = True
619  except KeyError:
620  pass
621 
622  if foundBkg:
623  bgStats = (bg[0].getStatsImage().getImage().array
624  for bg in bgList)
625  rec['skyBackground'] = sum(np.median(bg[np.isfinite(bg)]) for bg in bgStats)
626  else:
627  self.log.warning('Sky background not found for visit %d / ccd %d' %
628  (visit, det))
629  rec['skyBackground'] = -1.0
630  else:
631  rec['skyBackground'] = -1.0
632 
633  rec['used'] = 1
634 
635  def _makeSourceMapper(self, sourceSchema):
636  """
637  Make a schema mapper for fgcm sources
638 
639  Parameters
640  ----------
641  sourceSchema: `afwTable.Schema`
642  Default source schema from the butler
643 
644  Returns
645  -------
646  sourceMapper: `afwTable.schemaMapper`
647  Mapper to the FGCM source schema
648  """
649 
650  # create a mapper to the preferred output
651  sourceMapper = afwTable.SchemaMapper(sourceSchema)
652 
653  # map to ra/dec
654  sourceMapper.addMapping(sourceSchema['coord_ra'].asKey(), 'ra')
655  sourceMapper.addMapping(sourceSchema['coord_dec'].asKey(), 'dec')
656  sourceMapper.addMapping(sourceSchema['slot_Centroid_x'].asKey(), 'x')
657  sourceMapper.addMapping(sourceSchema['slot_Centroid_y'].asKey(), 'y')
658  # Add the mapping if the field exists in the input catalog.
659  # If the field does not exist, simply add it (set to False).
660  # This field is not required for calibration, but is useful
661  # to collate if available.
662  try:
663  sourceMapper.addMapping(sourceSchema[self.config.psfCandidateName].asKey(),
664  'psf_candidate')
665  except LookupError:
666  sourceMapper.editOutputSchema().addField(
667  "psf_candidate", type='Flag',
668  doc=("Flag set if the source was a candidate for PSF determination, "
669  "as determined by the star selector."))
670 
671  # and add the fields we want
672  sourceMapper.editOutputSchema().addField(
673  "visit", type=np.int32, doc="Visit number")
674  sourceMapper.editOutputSchema().addField(
675  "ccd", type=np.int32, doc="CCD number")
676  sourceMapper.editOutputSchema().addField(
677  "instMag", type=np.float32, doc="Instrumental magnitude")
678  sourceMapper.editOutputSchema().addField(
679  "instMagErr", type=np.float32, doc="Instrumental magnitude error")
680  sourceMapper.editOutputSchema().addField(
681  "jacobian", type=np.float32, doc="Relative pixel scale from wcs jacobian")
682  sourceMapper.editOutputSchema().addField(
683  "deltaMagBkg", type=np.float32, doc="Change in magnitude due to local background offset")
684 
685  return sourceMapper
686 
687  def fgcmMatchStars(self, visitCat, obsCat, lutDataRef=None):
688  """
689  Use FGCM code to match observations into unique stars.
690 
691  Parameters
692  ----------
693  visitCat: `afw.table.BaseCatalog`
694  Catalog with visit data for fgcm
695  obsCat: `afw.table.BaseCatalog`
696  Full catalog of star observations for fgcm
697  lutDataRef: `lsst.daf.persistence.ButlerDataRef` or
698  `lsst.daf.butler.DeferredDatasetHandle`, optional
699  Data reference to fgcm look-up table (used if matching reference stars).
700 
701  Returns
702  -------
703  fgcmStarIdCat: `afw.table.BaseCatalog`
704  Catalog of unique star identifiers and index keys
705  fgcmStarIndicesCat: `afwTable.BaseCatalog`
706  Catalog of unique star indices
707  fgcmRefCat: `afw.table.BaseCatalog`
708  Catalog of matched reference stars.
709  Will be None if `config.doReferenceMatches` is False.
710  """
711  # get filter names into a numpy array...
712  # This is the type that is expected by the fgcm code
713  visitFilterNames = np.zeros(len(visitCat), dtype='a30')
714  for i in range(len(visitCat)):
715  visitFilterNames[i] = visitCat[i]['physicalFilter']
716 
717  # match to put filterNames with observations
718  visitIndex = np.searchsorted(visitCat['visit'],
719  obsCat['visit'])
720 
721  obsFilterNames = visitFilterNames[visitIndex]
722 
723  if self.config.doReferenceMatches:
724  # Get the reference filter names, using the LUT
725  lutCat = lutDataRef.get()
726 
727  stdFilterDict = {filterName: stdFilter for (filterName, stdFilter) in
728  zip(lutCat[0]['physicalFilters'].split(','),
729  lutCat[0]['stdPhysicalFilters'].split(','))}
730  stdLambdaDict = {stdFilter: stdLambda for (stdFilter, stdLambda) in
731  zip(lutCat[0]['stdPhysicalFilters'].split(','),
732  lutCat[0]['lambdaStdFilter'])}
733 
734  del lutCat
735 
736  referenceFilterNames = self._getReferenceFilterNames_getReferenceFilterNames(visitCat,
737  stdFilterDict,
738  stdLambdaDict)
739  self.log.info("Using the following reference filters: %s" %
740  (', '.join(referenceFilterNames)))
741 
742  else:
743  # This should be an empty list
744  referenceFilterNames = []
745 
746  # make the fgcm starConfig dict
747  starConfig = {'logger': self.log,
748  'filterToBand': self.config.physicalFilterMap,
749  'requiredBands': self.config.requiredBands,
750  'minPerBand': self.config.minPerBand,
751  'matchRadius': self.config.matchRadius,
752  'isolationRadius': self.config.isolationRadius,
753  'matchNSide': self.config.matchNside,
754  'coarseNSide': self.config.coarseNside,
755  'densNSide': self.config.densityCutNside,
756  'densMaxPerPixel': self.config.densityCutMaxPerPixel,
757  'randomSeed': self.config.randomSeed,
758  'primaryBands': self.config.primaryBands,
759  'referenceFilterNames': referenceFilterNames}
760 
761  # initialize the FgcmMakeStars object
762  fgcmMakeStars = fgcm.FgcmMakeStars(starConfig)
763 
764  # make the primary stars
765  # note that the ra/dec native Angle format is radians
766  # We determine the conversion from the native units (typically
767  # radians) to degrees for the first observation. This allows us
768  # to treate ra/dec as numpy arrays rather than Angles, which would
769  # be approximately 600x slower.
770  conv = obsCat[0]['ra'].asDegrees() / float(obsCat[0]['ra'])
771  fgcmMakeStars.makePrimaryStars(obsCat['ra'] * conv,
772  obsCat['dec'] * conv,
773  filterNameArray=obsFilterNames,
774  bandSelected=False)
775 
776  # and match all the stars
777  fgcmMakeStars.makeMatchedStars(obsCat['ra'] * conv,
778  obsCat['dec'] * conv,
779  obsFilterNames)
780 
781  if self.config.doReferenceMatches:
782  fgcmMakeStars.makeReferenceMatches(self.fgcmLoadReferenceCatalog)
783 
784  # now persist
785 
786  objSchema = self._makeFgcmObjSchema_makeFgcmObjSchema()
787 
788  # make catalog and records
789  fgcmStarIdCat = afwTable.BaseCatalog(objSchema)
790  fgcmStarIdCat.reserve(fgcmMakeStars.objIndexCat.size)
791  for i in range(fgcmMakeStars.objIndexCat.size):
792  fgcmStarIdCat.addNew()
793 
794  # fill the catalog
795  fgcmStarIdCat['fgcm_id'][:] = fgcmMakeStars.objIndexCat['fgcm_id']
796  fgcmStarIdCat['ra'][:] = fgcmMakeStars.objIndexCat['ra']
797  fgcmStarIdCat['dec'][:] = fgcmMakeStars.objIndexCat['dec']
798  fgcmStarIdCat['obsArrIndex'][:] = fgcmMakeStars.objIndexCat['obsarrindex']
799  fgcmStarIdCat['nObs'][:] = fgcmMakeStars.objIndexCat['nobs']
800 
801  obsSchema = self._makeFgcmObsSchema_makeFgcmObsSchema()
802 
803  fgcmStarIndicesCat = afwTable.BaseCatalog(obsSchema)
804  fgcmStarIndicesCat.reserve(fgcmMakeStars.obsIndexCat.size)
805  for i in range(fgcmMakeStars.obsIndexCat.size):
806  fgcmStarIndicesCat.addNew()
807 
808  fgcmStarIndicesCat['obsIndex'][:] = fgcmMakeStars.obsIndexCat['obsindex']
809 
810  if self.config.doReferenceMatches:
811  refSchema = self._makeFgcmRefSchema_makeFgcmRefSchema(len(referenceFilterNames))
812 
813  fgcmRefCat = afwTable.BaseCatalog(refSchema)
814  fgcmRefCat.reserve(fgcmMakeStars.referenceCat.size)
815 
816  for i in range(fgcmMakeStars.referenceCat.size):
817  fgcmRefCat.addNew()
818 
819  fgcmRefCat['fgcm_id'][:] = fgcmMakeStars.referenceCat['fgcm_id']
820  fgcmRefCat['refMag'][:, :] = fgcmMakeStars.referenceCat['refMag']
821  fgcmRefCat['refMagErr'][:, :] = fgcmMakeStars.referenceCat['refMagErr']
822 
823  md = PropertyList()
824  md.set("REFSTARS_FORMAT_VERSION", REFSTARS_FORMAT_VERSION)
825  md.set("FILTERNAMES", referenceFilterNames)
826  fgcmRefCat.setMetadata(md)
827 
828  else:
829  fgcmRefCat = None
830 
831  return fgcmStarIdCat, fgcmStarIndicesCat, fgcmRefCat
832 
833  def _makeFgcmVisitSchema(self, nCcd):
834  """
835  Make a schema for an fgcmVisitCatalog
836 
837  Parameters
838  ----------
839  nCcd: `int`
840  Number of CCDs in the camera
841 
842  Returns
843  -------
844  schema: `afwTable.Schema`
845  """
846 
847  schema = afwTable.Schema()
848  schema.addField('visit', type=np.int32, doc="Visit number")
849  schema.addField('physicalFilter', type=str, size=30, doc="Physical filter")
850  schema.addField('telra', type=np.float64, doc="Pointing RA (deg)")
851  schema.addField('teldec', type=np.float64, doc="Pointing Dec (deg)")
852  schema.addField('telha', type=np.float64, doc="Pointing Hour Angle (deg)")
853  schema.addField('telrot', type=np.float64, doc="Camera rotation (deg)")
854  schema.addField('mjd', type=np.float64, doc="MJD of visit")
855  schema.addField('exptime', type=np.float32, doc="Exposure time")
856  schema.addField('pmb', type=np.float32, doc="Pressure (millibar)")
857  schema.addField('psfSigma', type=np.float32, doc="PSF sigma (reference CCD)")
858  schema.addField('deltaAper', type=np.float32, doc="Delta-aperture")
859  schema.addField('skyBackground', type=np.float32, doc="Sky background (ADU) (reference CCD)")
860  # the following field is not used yet
861  schema.addField('deepFlag', type=np.int32, doc="Deep observation")
862  schema.addField('scaling', type='ArrayD', doc="Scaling applied due to flat adjustment",
863  size=nCcd)
864  schema.addField('used', type=np.int32, doc="This visit has been ingested.")
865  schema.addField('sources_read', type='Flag', doc="This visit had sources read.")
866 
867  return schema
868 
869  def _makeFgcmObjSchema(self):
870  """
871  Make a schema for the objIndexCat from fgcmMakeStars
872 
873  Returns
874  -------
875  schema: `afwTable.Schema`
876  """
877 
878  objSchema = afwTable.Schema()
879  objSchema.addField('fgcm_id', type=np.int32, doc='FGCM Unique ID')
880  # Will investigate making these angles...
881  objSchema.addField('ra', type=np.float64, doc='Mean object RA (deg)')
882  objSchema.addField('dec', type=np.float64, doc='Mean object Dec (deg)')
883  objSchema.addField('obsArrIndex', type=np.int32,
884  doc='Index in obsIndexTable for first observation')
885  objSchema.addField('nObs', type=np.int32, doc='Total number of observations')
886 
887  return objSchema
888 
889  def _makeFgcmObsSchema(self):
890  """
891  Make a schema for the obsIndexCat from fgcmMakeStars
892 
893  Returns
894  -------
895  schema: `afwTable.Schema`
896  """
897 
898  obsSchema = afwTable.Schema()
899  obsSchema.addField('obsIndex', type=np.int32, doc='Index in observation table')
900 
901  return obsSchema
902 
903  def _makeFgcmRefSchema(self, nReferenceBands):
904  """
905  Make a schema for the referenceCat from fgcmMakeStars
906 
907  Parameters
908  ----------
909  nReferenceBands: `int`
910  Number of reference bands
911 
912  Returns
913  -------
914  schema: `afwTable.Schema`
915  """
916 
917  refSchema = afwTable.Schema()
918  refSchema.addField('fgcm_id', type=np.int32, doc='FGCM Unique ID')
919  refSchema.addField('refMag', type='ArrayF', doc='Reference magnitude array (AB)',
920  size=nReferenceBands)
921  refSchema.addField('refMagErr', type='ArrayF', doc='Reference magnitude error array',
922  size=nReferenceBands)
923 
924  return refSchema
925 
926  def _getReferenceFilterNames(self, visitCat, stdFilterDict, stdLambdaDict):
927  """
928  Get the reference filter names, in wavelength order, from the visitCat and
929  information from the look-up-table.
930 
931  Parameters
932  ----------
933  visitCat: `afw.table.BaseCatalog`
934  Catalog with visit data for FGCM
935  stdFilterDict: `dict`
936  Mapping of filterName to stdFilterName from LUT
937  stdLambdaDict: `dict`
938  Mapping of stdFilterName to stdLambda from LUT
939 
940  Returns
941  -------
942  referenceFilterNames: `list`
943  Wavelength-ordered list of reference filter names
944  """
945 
946  # Find the unique list of filter names in visitCat
947  filterNames = np.unique(visitCat.asAstropy()['physicalFilter'])
948 
949  # Find the unique list of "standard" filters
950  stdFilterNames = {stdFilterDict[filterName] for filterName in filterNames}
951 
952  # And sort these by wavelength
953  referenceFilterNames = sorted(stdFilterNames, key=stdLambdaDict.get)
954 
955  return referenceFilterNames
Defines the fields and offsets for a table.
Definition: Schema.h:51
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
def fgcmMatchStars(self, visitCat, obsCat, lutDataRef=None)
def _getReferenceFilterNames(self, visitCat, stdFilterDict, stdLambdaDict)
def _fillVisitCatalog(self, visitCat, groupedDataRefs, bkgDataRefDict=None, visitCatDataRef=None)
def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat, sourceSchema, camera, calibFluxApertureRadius=None, visitCatDataRef=None, starObsDataRef=None, inStarObsCat=None)
def __init__(self, initInputs=None, butler=None, **kwargs)
def fgcmMakeVisitCatalog(self, camera, groupedDataRefs, bkgDataRefDict=None, visitCatDataRef=None, inVisitCat=None)
def _findAndGroupDataRefsGen2(self, butler, camera, dataRefs)
An integer coordinate rectangle.
Definition: Box.h:55
daf::base::PropertyList * list
Definition: fits.cc:913
def computeApertureRadiusFromDataRef(dataRef, fluxField)
Definition: utilities.py:812