LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
forcedPhotCcd.py
Go to the documentation of this file.
1 # This file is part of meas_base.
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 
22 import collections
23 import logging
24 import pandas as pd
25 import numpy as np
26 
27 import lsst.pex.config
29 import lsst.pipe.base
30 import lsst.geom
31 import lsst.afw.detection
32 import lsst.afw.geom
33 import lsst.afw.image
34 import lsst.afw.table
35 import lsst.sphgeom
36 
37 from lsst.obs.base import ExposureIdInfo
38 from lsst.pipe.base import PipelineTaskConnections
39 import lsst.pipe.base.connectionTypes as cT
40 
41 import lsst.pipe.base as pipeBase
42 from lsst.skymap import BaseSkyMap
43 
44 from .references import MultiBandReferencesTask
45 from .forcedMeasurement import ForcedMeasurementTask
46 from .applyApCorr import ApplyApCorrTask
47 from .catalogCalculation import CatalogCalculationTask
48 
49 try:
50  from lsst.meas.mosaic import applyMosaicResults
51 except ImportError:
52  applyMosaicResults = None
53 
54 __all__ = ("PerTractCcdDataIdContainer", "ForcedPhotCcdConfig", "ForcedPhotCcdTask", "imageOverlapsTract",
55  "ForcedPhotCcdFromDataFrameTask", "ForcedPhotCcdFromDataFrameConfig")
56 
57 
58 class PerTractCcdDataIdContainer(pipeBase.DataIdContainer):
59  """A data ID container which combines raw data IDs with a tract.
60 
61  Notes
62  -----
63  Required because we need to add "tract" to the raw data ID keys (defined as
64  whatever we use for ``src``) when no tract is provided (so that the user is
65  not required to know which tracts are spanned by the raw data ID).
66 
67  This subclass of `~lsst.pipe.base.DataIdContainer` assumes that a calexp is
68  being measured using the detection information, a set of reference
69  catalogs, from the set of coadds which intersect with the calexp. It needs
70  the calexp id (e.g. visit, raft, sensor), but is also uses the tract to
71  decide what set of coadds to use. The references from the tract whose
72  patches intersect with the calexp are used.
73  """
74 
75  def makeDataRefList(self, namespace):
76  """Make self.refList from self.idList
77  """
78  if self.datasetType is None:
79  raise RuntimeError("Must call setDatasetType first")
80  log = logging.getLogger("meas.base.forcedPhotCcd.PerTractCcdDataIdContainer")
81  skymap = None
82  visitTract = collections.defaultdict(set) # Set of tracts for each visit
83  visitRefs = collections.defaultdict(list) # List of data references for each visit
84  for dataId in self.idList:
85  if "tract" not in dataId:
86  # Discover which tracts the data overlaps
87  log.info("Reading WCS for components of dataId=%s to determine tracts", dict(dataId))
88  if skymap is None:
89  skymap = namespace.butler.get(namespace.config.coaddName + "Coadd_skyMap")
90 
91  for ref in namespace.butler.subset("calexp", dataId=dataId):
92  if not ref.datasetExists("calexp"):
93  continue
94 
95  visit = ref.dataId["visit"]
96  visitRefs[visit].append(ref)
97 
98  md = ref.get("calexp_md", immediate=True)
99  wcs = lsst.afw.geom.makeSkyWcs(md)
101  # Going with just the nearest tract. Since we're throwing all tracts for the visit
102  # together, this shouldn't be a problem unless the tracts are much smaller than a CCD.
103  tract = skymap.findTract(wcs.pixelToSky(box.getCenter()))
104  if imageOverlapsTract(tract, wcs, box):
105  visitTract[visit].add(tract.getId())
106  else:
107  self.refList.extend(ref for ref in namespace.butler.subset(self.datasetType, dataId=dataId))
108 
109  # Ensure all components of a visit are kept together by putting them all in the same set of tracts
110  for visit, tractSet in visitTract.items():
111  for ref in visitRefs[visit]:
112  for tract in tractSet:
113  self.refList.append(namespace.butler.dataRef(datasetType=self.datasetType,
114  dataId=ref.dataId, tract=tract))
115  if visitTract:
116  tractCounter = collections.Counter()
117  for tractSet in visitTract.values():
118  tractCounter.update(tractSet)
119  log.info("Number of visits for each tract: %s", dict(tractCounter))
120 
121 
122 def imageOverlapsTract(tract, imageWcs, imageBox):
123  """Return whether the given bounding box overlaps the tract given a WCS.
124 
125  Parameters
126  ----------
127  tract : `lsst.skymap.TractInfo`
128  TractInfo specifying a tract.
129  imageWcs : `lsst.afw.geom.SkyWcs`
130  World coordinate system for the image.
131  imageBox : `lsst.geom.Box2I`
132  Bounding box for the image.
133 
134  Returns
135  -------
136  overlap : `bool`
137  `True` if the bounding box overlaps the tract; `False` otherwise.
138  """
139  tractPoly = tract.getOuterSkyPolygon()
140 
141  imagePixelCorners = lsst.geom.Box2D(imageBox).getCorners()
142  try:
143  imageSkyCorners = imageWcs.pixelToSky(imagePixelCorners)
144  except lsst.pex.exceptions.LsstCppException as e:
145  # Protecting ourselves from awful Wcs solutions in input images
146  if (not isinstance(e.message, lsst.pex.exceptions.DomainErrorException)
147  and not isinstance(e.message, lsst.pex.exceptions.RuntimeErrorException)):
148  raise
149  return False
150 
151  imagePoly = lsst.sphgeom.ConvexPolygon.convexHull([coord.getVector() for coord in imageSkyCorners])
152  return tractPoly.intersects(imagePoly) # "intersects" also covers "contains" or "is contained by"
153 
154 
155 class ForcedPhotCcdConnections(PipelineTaskConnections,
156  dimensions=("instrument", "visit", "detector", "skymap", "tract"),
157  defaultTemplates={"inputCoaddName": "deep",
158  "inputName": "calexp"}):
159  inputSchema = cT.InitInput(
160  doc="Schema for the input measurement catalogs.",
161  name="{inputCoaddName}Coadd_ref_schema",
162  storageClass="SourceCatalog",
163  )
164  outputSchema = cT.InitOutput(
165  doc="Schema for the output forced measurement catalogs.",
166  name="forced_src_schema",
167  storageClass="SourceCatalog",
168  )
169  exposure = cT.Input(
170  doc="Input exposure to perform photometry on.",
171  name="{inputName}",
172  storageClass="ExposureF",
173  dimensions=["instrument", "visit", "detector"],
174  )
175  refCat = cT.Input(
176  doc="Catalog of shapes and positions at which to force photometry.",
177  name="{inputCoaddName}Coadd_ref",
178  storageClass="SourceCatalog",
179  dimensions=["skymap", "tract", "patch"],
180  multiple=True,
181  deferLoad=True,
182  )
183  skyMap = cT.Input(
184  doc="SkyMap dataset that defines the coordinate system of the reference catalog.",
185  name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
186  storageClass="SkyMap",
187  dimensions=["skymap"],
188  )
189  measCat = cT.Output(
190  doc="Output forced photometry catalog.",
191  name="forced_src",
192  storageClass="SourceCatalog",
193  dimensions=["instrument", "visit", "detector", "skymap", "tract"],
194  )
195 
196 
197 class ForcedPhotCcdConfig(pipeBase.PipelineTaskConfig,
198  pipelineConnections=ForcedPhotCcdConnections):
199  """Config class for forced measurement driver task."""
201  target=MultiBandReferencesTask,
202  doc="subtask to retrieve reference source catalog"
203  )
204  measurement = lsst.pex.config.ConfigurableField(
205  target=ForcedMeasurementTask,
206  doc="subtask to do forced measurement"
207  )
208  coaddName = lsst.pex.config.Field(
209  doc="coadd name: typically one of deep or goodSeeing",
210  dtype=str,
211  default="deep",
212  )
213  doApCorr = lsst.pex.config.Field(
214  dtype=bool,
215  default=True,
216  doc="Run subtask to apply aperture corrections"
217  )
218  applyApCorr = lsst.pex.config.ConfigurableField(
219  target=ApplyApCorrTask,
220  doc="Subtask to apply aperture corrections"
221  )
222  catalogCalculation = lsst.pex.config.ConfigurableField(
223  target=CatalogCalculationTask,
224  doc="Subtask to run catalogCalculation plugins on catalog"
225  )
226  doApplyUberCal = lsst.pex.config.Field(
227  dtype=bool,
228  doc="Apply meas_mosaic ubercal results to input calexps?",
229  default=False,
230  deprecated="Deprecated by DM-23352; use doApplyExternalPhotoCalib and doApplyExternalSkyWcs instead",
231  )
232  doApplyExternalPhotoCalib = lsst.pex.config.Field(
233  dtype=bool,
234  default=False,
235  doc=("Whether to apply external photometric calibration via an "
236  "`lsst.afw.image.PhotoCalib` object. Uses the "
237  "``externalPhotoCalibName`` field to determine which calibration "
238  "to load."),
239  )
240  doApplyExternalSkyWcs = lsst.pex.config.Field(
241  dtype=bool,
242  default=False,
243  doc=("Whether to apply external astrometric calibration via an "
244  "`lsst.afw.geom.SkyWcs` object. Uses ``externalSkyWcsName`` "
245  "field to determine which calibration to load."),
246  )
247  doApplySkyCorr = lsst.pex.config.Field(
248  dtype=bool,
249  default=False,
250  doc="Apply sky correction?",
251  )
252  includePhotoCalibVar = lsst.pex.config.Field(
253  dtype=bool,
254  default=False,
255  doc="Add photometric calibration variance to warp variance plane?",
256  )
257  externalPhotoCalibName = lsst.pex.config.ChoiceField(
258  dtype=str,
259  doc=("Type of external PhotoCalib if ``doApplyExternalPhotoCalib`` is True. "
260  "Unused for Gen3 middleware."),
261  default="jointcal",
262  allowed={
263  "jointcal": "Use jointcal_photoCalib",
264  "fgcm": "Use fgcm_photoCalib",
265  "fgcm_tract": "Use fgcm_tract_photoCalib"
266  },
267  )
268  externalSkyWcsName = lsst.pex.config.ChoiceField(
269  dtype=str,
270  doc="Type of external SkyWcs if ``doApplyExternalSkyWcs`` is True. Unused for Gen3 middleware.",
271  default="jointcal",
272  allowed={
273  "jointcal": "Use jointcal_wcs"
274  },
275  )
276  footprintSource = lsst.pex.config.ChoiceField(
277  dtype=str,
278  doc="Where to obtain footprints to install in the measurement catalog, prior to measurement.",
279  allowed={
280  "transformed": "Transform footprints from the reference catalog (downgrades HeavyFootprints).",
281  "psf": ("Use the scaled shape of the PSF at the position of each source (does not generate "
282  "HeavyFootprints)."),
283  },
284  optional=True,
285  default="transformed",
286  )
287  psfFootprintScaling = lsst.pex.config.Field(
288  dtype=float,
289  doc="Scaling factor to apply to the PSF shape when footprintSource='psf' (ignored otherwise).",
290  default=3.0,
291  )
292 
293  def setDefaults(self):
294  # Docstring inherited.
295  # Make catalogCalculation a no-op by default as no modelFlux is setup by default in
296  # ForcedMeasurementTask
297  super().setDefaults()
298  self.measurement.plugins.names |= ['base_LocalPhotoCalib', 'base_LocalWcs']
299  self.catalogCalculation.plugins.names = []
300 
301 
302 class ForcedPhotCcdTask(pipeBase.PipelineTask, pipeBase.CmdLineTask):
303  """A command-line driver for performing forced measurement on CCD images.
304 
305  Parameters
306  ----------
307  butler : `lsst.daf.persistence.butler.Butler`, optional
308  A Butler which will be passed to the references subtask to allow it to
309  load its schema from disk. Optional, but must be specified if
310  ``refSchema`` is not; if both are specified, ``refSchema`` takes
311  precedence.
312  refSchema : `lsst.afw.table.Schema`, optional
313  The schema of the reference catalog, passed to the constructor of the
314  references subtask. Optional, but must be specified if ``butler`` is
315  not; if both are specified, ``refSchema`` takes precedence.
316  **kwds
317  Keyword arguments are passed to the supertask constructor.
318 
319  Notes
320  -----
321  The `runDataRef` method takes a `~lsst.daf.persistence.ButlerDataRef` argument
322  that corresponds to a single CCD. This should contain the data ID keys that
323  correspond to the ``forced_src`` dataset (the output dataset for this
324  task), which are typically all those used to specify the ``calexp`` dataset
325  (``visit``, ``raft``, ``sensor`` for LSST data) as well as a coadd tract.
326  The tract is used to look up the appropriate coadd measurement catalogs to
327  use as references (e.g. ``deepCoadd_src``; see
328  :lsst-task:`lsst.meas.base.references.CoaddSrcReferencesTask` for more
329  information). While the tract must be given as part of the dataRef, the
330  patches are determined automatically from the bounding box and WCS of the
331  calexp to be measured, and the filter used to fetch references is set via
332  the ``filter`` option in the configuration of
333  :lsst-task:`lsst.meas.base.references.BaseReferencesTask`).
334  """
335 
336  ConfigClass = ForcedPhotCcdConfig
337  RunnerClass = pipeBase.ButlerInitializedTaskRunner
338  _DefaultName = "forcedPhotCcd"
339  dataPrefix = ""
340 
341  def __init__(self, butler=None, refSchema=None, initInputs=None, **kwds):
342  super().__init__(**kwds)
343 
344  if initInputs is not None:
345  refSchema = initInputs['inputSchema'].schema
346 
347  self.makeSubtask("references", butler=butler, schema=refSchema)
348  if refSchema is None:
349  refSchema = self.references.schema
350  self.makeSubtask("measurement", refSchema=refSchema)
351  # It is necessary to get the schema internal to the forced measurement task until such a time
352  # that the schema is not owned by the measurement task, but is passed in by an external caller
353  if self.config.doApCorr:
354  self.makeSubtask("applyApCorr", schema=self.measurement.schema)
355  self.makeSubtask('catalogCalculation', schema=self.measurement.schema)
356  self.outputSchema = lsst.afw.table.SourceCatalog(self.measurement.schema)
357 
358  def runQuantum(self, butlerQC, inputRefs, outputRefs):
359  inputs = butlerQC.get(inputRefs)
360 
361  tract = butlerQC.quantum.dataId['tract']
362  skyMap = inputs.pop("skyMap")
363  inputs['refWcs'] = skyMap[tract].getWcs()
364 
365  inputs['refCat'] = self.mergeAndFilterReferences(inputs['exposure'], inputs['refCat'],
366  inputs['refWcs'])
367 
368  inputs['measCat'], inputs['exposureId'] = self.generateMeasCat(inputRefs.exposure.dataId,
369  inputs['exposure'],
370  inputs['refCat'], inputs['refWcs'],
371  "visit_detector")
372  self.attachFootprints(inputs['measCat'], inputs['refCat'], inputs['exposure'], inputs['refWcs'])
373  # TODO: apply external calibrations (DM-17062)
374  outputs = self.run(**inputs)
375  butlerQC.put(outputs, outputRefs)
376 
377  def mergeAndFilterReferences(self, exposure, refCats, refWcs):
378  """Filter reference catalog so that all sources are within the
379  boundaries of the exposure.
380 
381  Parameters
382  ----------
383  exposure : `lsst.afw.image.exposure.Exposure`
384  Exposure to generate the catalog for.
385  refCats : sequence of `lsst.daf.butler.DeferredDatasetHandle`
386  Handles for catalogs of shapes and positions at which to force
387  photometry.
388  refWcs : `lsst.afw.image.SkyWcs`
389  Reference world coordinate system.
390 
391  Returns
392  -------
393  refSources : `lsst.afw.table.SourceCatalog`
394  Filtered catalog of forced sources to measure.
395 
396  Notes
397  -----
398  Filtering the reference catalog is currently handled by Gen2
399  specific methods. To function for Gen3, this method copies
400  code segments to do the filtering and transformation. The
401  majority of this code is based on the methods of
402  lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader
403 
404  """
405 
406  # Step 1: Determine bounds of the exposure photometry will
407  # be performed on.
408  expWcs = exposure.getWcs()
409  expRegion = exposure.getBBox(lsst.afw.image.PARENT)
410  expBBox = lsst.geom.Box2D(expRegion)
411  expBoxCorners = expBBox.getCorners()
412  expSkyCorners = [expWcs.pixelToSky(corner).getVector() for
413  corner in expBoxCorners]
414  expPolygon = lsst.sphgeom.ConvexPolygon(expSkyCorners)
415 
416  # Step 2: Filter out reference catalog sources that are
417  # not contained within the exposure boundaries, or whose
418  # parents are not within the exposure boundaries. Note
419  # that within a single input refCat, the parents always
420  # appear before the children.
421  mergedRefCat = None
422  for refCat in refCats:
423  refCat = refCat.get()
424  if mergedRefCat is None:
425  mergedRefCat = lsst.afw.table.SourceCatalog(refCat.table)
426  containedIds = {0} # zero as a parent ID means "this is a parent"
427  for record in refCat:
428  if expPolygon.contains(record.getCoord().getVector()) and record.getParent() in containedIds:
429  record.setFootprint(record.getFootprint())
430  mergedRefCat.append(record)
431  containedIds.add(record.getId())
432  if mergedRefCat is None:
433  raise RuntimeError("No reference objects for forced photometry.")
434  mergedRefCat.sort(lsst.afw.table.SourceTable.getParentKey())
435  return mergedRefCat
436 
437  def generateMeasCat(self, exposureDataId, exposure, refCat, refWcs, idPackerName):
438  """Generate a measurement catalog for Gen3.
439 
440  Parameters
441  ----------
442  exposureDataId : `DataId`
443  Butler dataId for this exposure.
444  exposure : `lsst.afw.image.exposure.Exposure`
445  Exposure to generate the catalog for.
446  refCat : `lsst.afw.table.SourceCatalog`
447  Catalog of shapes and positions at which to force photometry.
448  refWcs : `lsst.afw.image.SkyWcs`
449  Reference world coordinate system.
450  idPackerName : `str`
451  Type of ID packer to construct from the registry.
452 
453  Returns
454  -------
455  measCat : `lsst.afw.table.SourceCatalog`
456  Catalog of forced sources to measure.
457  expId : `int`
458  Unique binary id associated with the input exposure
459  """
460  exposureIdInfo = ExposureIdInfo.fromDataId(exposureDataId, idPackerName)
461  idFactory = exposureIdInfo.makeSourceIdFactory()
462 
463  measCat = self.measurement.generateMeasCat(exposure, refCat, refWcs,
464  idFactory=idFactory)
465  return measCat, exposureIdInfo.expId
466 
467  def runDataRef(self, dataRef, psfCache=None):
468  """Perform forced measurement on a single exposure.
469 
470  Parameters
471  ----------
472  dataRef : `lsst.daf.persistence.ButlerDataRef`
473  Passed to the ``references`` subtask to obtain the reference WCS,
474  the ``getExposure`` method (implemented by derived classes) to
475  read the measurment image, and the ``fetchReferences`` method to
476  get the exposure and load the reference catalog (see
477  :lsst-task`lsst.meas.base.references.CoaddSrcReferencesTask`).
478  Refer to derived class documentation for details of the datasets
479  and data ID keys which are used.
480  psfCache : `int`, optional
481  Size of PSF cache, or `None`. The size of the PSF cache can have
482  a significant effect upon the runtime for complicated PSF models.
483 
484  Notes
485  -----
486  Sources are generated with ``generateMeasCat`` in the ``measurement``
487  subtask. These are passed to ``measurement``'s ``run`` method, which
488  fills the source catalog with the forced measurement results. The
489  sources are then passed to the ``writeOutputs`` method (implemented by
490  derived classes) which writes the outputs.
491  """
492  refWcs = self.references.getWcs(dataRef)
493  exposure = self.getExposure(dataRef)
494  if psfCache is not None:
495  exposure.getPsf().setCacheSize(psfCache)
496  refCat = self.fetchReferences(dataRef, exposure)
497  measCat = self.measurement.generateMeasCat(exposure, refCat, refWcs,
498  idFactory=self.makeIdFactory(dataRef))
499  self.log.info("Performing forced measurement on %s", dataRef.dataId)
500  self.attachFootprints(measCat, refCat, exposure, refWcs)
501 
502  exposureId = self.getExposureId(dataRef)
503 
504  forcedPhotResult = self.run(measCat, exposure, refCat, refWcs, exposureId=exposureId)
505 
506  self.writeOutput(dataRef, forcedPhotResult.measCat)
507 
508  def run(self, measCat, exposure, refCat, refWcs, exposureId=None):
509  """Perform forced measurement on a single exposure.
510 
511  Parameters
512  ----------
513  measCat : `lsst.afw.table.SourceCatalog`
514  The measurement catalog, based on the sources listed in the
515  reference catalog.
516  exposure : `lsst.afw.image.Exposure`
517  The measurement image upon which to perform forced detection.
518  refCat : `lsst.afw.table.SourceCatalog`
519  The reference catalog of sources to measure.
520  refWcs : `lsst.afw.image.SkyWcs`
521  The WCS for the references.
522  exposureId : `int`
523  Optional unique exposureId used for random seed in measurement
524  task.
525 
526  Returns
527  -------
528  result : `lsst.pipe.base.Struct`
529  Structure with fields:
530 
531  ``measCat``
532  Catalog of forced measurement results
533  (`lsst.afw.table.SourceCatalog`).
534  """
535  self.measurement.run(measCat, exposure, refCat, refWcs, exposureId=exposureId)
536  if self.config.doApCorr:
537  self.applyApCorr.run(
538  catalog=measCat,
539  apCorrMap=exposure.getInfo().getApCorrMap()
540  )
541  self.catalogCalculation.run(measCat)
542 
543  return pipeBase.Struct(measCat=measCat)
544 
545  def makeIdFactory(self, dataRef):
546  """Create an object that generates globally unique source IDs.
547 
548  Source IDs are created based on a per-CCD ID and the ID of the CCD
549  itself.
550 
551  Parameters
552  ----------
553  dataRef : `lsst.daf.persistence.ButlerDataRef`
554  Butler data reference. The ``ccdExposureId_bits`` and
555  ``ccdExposureId`` datasets are accessed. The data ID must have the
556  keys that correspond to ``ccdExposureId``, which are generally the
557  same as those that correspond to ``calexp`` (``visit``, ``raft``,
558  ``sensor`` for LSST data).
559  """
560  exposureIdInfo = ExposureIdInfo(int(dataRef.get("ccdExposureId")), dataRef.get("ccdExposureId_bits"))
561  return exposureIdInfo.makeSourceIdFactory()
562 
563  def getExposureId(self, dataRef):
564  return int(dataRef.get("ccdExposureId", immediate=True))
565 
566  def fetchReferences(self, dataRef, exposure):
567  """Get sources that overlap the exposure.
568 
569  Parameters
570  ----------
571  dataRef : `lsst.daf.persistence.ButlerDataRef`
572  Butler data reference corresponding to the image to be measured;
573  should have ``tract``, ``patch``, and ``filter`` keys.
574  exposure : `lsst.afw.image.Exposure`
575  The image to be measured (used only to obtain a WCS and bounding
576  box).
577 
578  Returns
579  -------
580  referencs : `lsst.afw.table.SourceCatalog`
581  Catalog of sources that overlap the exposure
582 
583  Notes
584  -----
585  The returned catalog is sorted by ID and guarantees that all included
586  children have their parent included and that all Footprints are valid.
587 
588  All work is delegated to the references subtask; see
589  :lsst-task:`lsst.meas.base.references.CoaddSrcReferencesTask`
590  for information about the default behavior.
591  """
592  references = lsst.afw.table.SourceCatalog(self.references.schema)
593  badParents = set()
594  unfiltered = self.references.fetchInBox(dataRef, exposure.getBBox(), exposure.getWcs())
595  for record in unfiltered:
596  if record.getFootprint() is None or record.getFootprint().getArea() == 0:
597  if record.getParent() != 0:
598  self.log.warning("Skipping reference %s (child of %s) with bad Footprint",
599  record.getId(), record.getParent())
600  else:
601  self.log.warning("Skipping reference parent %s with bad Footprint", record.getId())
602  badParents.add(record.getId())
603  elif record.getParent() not in badParents:
604  references.append(record)
605  # catalog must be sorted by parent ID for lsst.afw.table.getChildren to work
606  references.sort(lsst.afw.table.SourceTable.getParentKey())
607  return references
608 
609  def attachFootprints(self, sources, refCat, exposure, refWcs):
610  r"""Attach footprints to blank sources prior to measurements.
611 
612  Notes
613  -----
614  `~lsst.afw.detection.Footprint`\ s for forced photometry must be in the
615  pixel coordinate system of the image being measured, while the actual
616  detections may start out in a different coordinate system.
617 
618  Subclasses of this class may implement this method to define how
619  those `~lsst.afw.detection.Footprint`\ s should be generated.
620 
621  This default implementation transforms depends on the
622  ``footprintSource`` configuration parameter.
623  """
624  if self.config.footprintSource == "transformed":
625  return self.measurement.attachTransformedFootprints(sources, refCat, exposure, refWcs)
626  elif self.config.footprintSource == "psf":
627  return self.measurement.attachPsfShapeFootprints(sources, exposure,
628  scaling=self.config.psfFootprintScaling)
629 
630  def getExposure(self, dataRef):
631  """Read input exposure for measurement.
632 
633  Parameters
634  ----------
635  dataRef : `lsst.daf.persistence.ButlerDataRef`
636  Butler data reference.
637  """
638  exposure = dataRef.get(self.dataPrefix + "calexp", immediate=True)
639 
640  if self.config.doApplyExternalPhotoCalib:
641  source = f"{self.config.externalPhotoCalibName}_photoCalib"
642  self.log.info("Applying external photoCalib from %s", source)
643  photoCalib = dataRef.get(source)
644  exposure.setPhotoCalib(photoCalib) # No need for calibrateImage; having the photoCalib suffices
645 
646  if self.config.doApplyExternalSkyWcs:
647  source = f"{self.config.externalSkyWcsName}_wcs"
648  self.log.info("Applying external skyWcs from %s", source)
649  skyWcs = dataRef.get(source)
650  exposure.setWcs(skyWcs)
651 
652  if self.config.doApplySkyCorr:
653  self.log.info("Apply sky correction")
654  skyCorr = dataRef.get("skyCorr")
655  exposure.maskedImage -= skyCorr.getImage()
656 
657  return exposure
658 
659  def writeOutput(self, dataRef, sources):
660  """Write forced source table
661 
662  Parameters
663  ----------
664  dataRef : `lsst.daf.persistence.ButlerDataRef`
665  Butler data reference. The forced_src dataset (with
666  self.dataPrefix prepended) is all that will be modified.
667  sources : `lsst.afw.table.SourceCatalog`
668  Catalog of sources to save.
669  """
670  dataRef.put(sources, self.dataPrefix + "forced_src", flags=lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS)
671 
672  def getSchemaCatalogs(self):
673  """The schema catalogs that will be used by this task.
674 
675  Returns
676  -------
677  schemaCatalogs : `dict`
678  Dictionary mapping dataset type to schema catalog.
679 
680  Notes
681  -----
682  There is only one schema for each type of forced measurement. The
683  dataset type for this measurement is defined in the mapper.
684  """
685  catalog = lsst.afw.table.SourceCatalog(self.measurement.schema)
686  catalog.getTable().setMetadata(self.measurement.algMetadata)
687  datasetType = self.dataPrefix + "forced_src"
688  return {datasetType: catalog}
689 
690  def _getConfigName(self):
691  # Documented in superclass.
692  return self.dataPrefix + "forcedPhotCcd_config"
693 
694  def _getMetadataName(self):
695  # Documented in superclass
696  return self.dataPrefix + "forcedPhotCcd_metadata"
697 
698  @classmethod
699  def _makeArgumentParser(cls):
700  parser = pipeBase.ArgumentParser(name=cls._DefaultName)
701  parser.add_id_argument("--id", "forced_src", help="data ID with raw CCD keys [+ tract optionally], "
702  "e.g. --id visit=12345 ccd=1,2 [tract=0]",
703  ContainerClass=PerTractCcdDataIdContainer)
704  return parser
705 
706 
707 class ForcedPhotCcdFromDataFrameConnections(PipelineTaskConnections,
708  dimensions=("instrument", "visit", "detector", "skymap", "tract"),
709  defaultTemplates={"inputCoaddName": "goodSeeing",
710  "inputName": "calexp"}):
711  refCat = cT.Input(
712  doc="Catalog of positions at which to force photometry.",
713  name="{inputCoaddName}Diff_fullDiaObjTable",
714  storageClass="DataFrame",
715  dimensions=["skymap", "tract", "patch"],
716  multiple=True,
717  deferLoad=True,
718  )
719  exposure = cT.Input(
720  doc="Input exposure to perform photometry on.",
721  name="{inputName}",
722  storageClass="ExposureF",
723  dimensions=["instrument", "visit", "detector"],
724  )
725  measCat = cT.Output(
726  doc="Output forced photometry catalog.",
727  name="forced_src_diaObject",
728  storageClass="SourceCatalog",
729  dimensions=["instrument", "visit", "detector", "skymap", "tract"],
730  )
731  outputSchema = cT.InitOutput(
732  doc="Schema for the output forced measurement catalogs.",
733  name="forced_src_diaObject_schema",
734  storageClass="SourceCatalog",
735  )
736 
737 
738 class ForcedPhotCcdFromDataFrameConfig(ForcedPhotCcdConfig,
739  pipelineConnections=ForcedPhotCcdFromDataFrameConnections):
740  def setDefaults(self):
741  super().setDefaults()
742  self.footprintSource = "psf"
743  self.measurement.doReplaceWithNoise = False
744  self.measurement.plugins = ["base_TransformedCentroidFromCoord", "base_PsfFlux", "base_PixelFlags"]
745  self.measurement.copyColumns = {'id': 'diaObjectId', 'coord_ra': 'coord_ra', 'coord_dec': 'coord_dec'}
746  self.measurement.slots.centroid = "base_TransformedCentroidFromCoord"
747  self.measurement.slots.psfFlux = "base_PsfFlux"
748  self.measurement.slots.shape = None
749 
750  def validate(self):
751  super().validate()
752  if self.footprintSource == "transformed":
753  raise ValueError("Cannot transform footprints from reference catalog, "
754  "because DataFrames can't hold footprints.")
755 
756 
757 class ForcedPhotCcdFromDataFrameTask(ForcedPhotCcdTask):
758  """Force Photometry on a per-detector exposure with coords from a DataFrame
759 
760  Uses input from a DataFrame instead of SourceCatalog
761  like the base class ForcedPhotCcd does.
762  Writes out a SourceCatalog so that the downstream
763  WriteForcedSourceTableTask can be reused with output from this Task.
764  """
765  _DefaultName = "forcedPhotCcdFromDataFrame"
766  ConfigClass = ForcedPhotCcdFromDataFrameConfig
767 
768  def __init__(self, butler=None, refSchema=None, initInputs=None, **kwds):
769  # Parent's init assumes that we have a reference schema; Cannot reuse
770  pipeBase.PipelineTask.__init__(self, **kwds)
771 
772  self.makeSubtask("measurement", refSchema=lsst.afw.table.SourceTable.makeMinimalSchema())
773 
774  if self.config.doApCorr:
775  self.makeSubtask("applyApCorr", schema=self.measurement.schema)
776  self.makeSubtask('catalogCalculation', schema=self.measurement.schema)
777  self.outputSchema = lsst.afw.table.SourceCatalog(self.measurement.schema)
778 
779  def runQuantum(self, butlerQC, inputRefs, outputRefs):
780  inputs = butlerQC.get(inputRefs)
781  self.log.info("Filtering ref cats: %s", ','.join([str(i.dataId) for i in inputs['refCat']]))
782  refCat = self.df2RefCat([i.get(parameters={"columns": ['diaObjectId', 'ra', 'decl']})
783  for i in inputs['refCat']],
784  inputs['exposure'].getBBox(), inputs['exposure'].getWcs())
785  inputs['refCat'] = refCat
786  inputs['refWcs'] = inputs['exposure'].getWcs()
787  inputs['measCat'], inputs['exposureId'] = self.generateMeasCat(inputRefs.exposure.dataId,
788  inputs['exposure'], inputs['refCat'],
789  inputs['refWcs'],
790  "visit_detector")
791  self.attachFootprints(inputs["measCat"], inputs["refCat"], inputs["exposure"], inputs["refWcs"])
792  outputs = self.run(**inputs)
793  butlerQC.put(outputs, outputRefs)
794 
795  def df2RefCat(self, dfList, exposureBBox, exposureWcs):
796  """Convert list of DataFrames to reference catalog
797 
798  Concatenate list of DataFrames presumably from multiple patches and
799  downselect rows that overlap the exposureBBox using the exposureWcs.
800 
801  Parameters
802  ----------
803  dfList : `list` of `pandas.DataFrame`
804  Each element containst diaObjects with ra/decl position in degrees
805  Columns 'diaObjectId', 'ra', 'decl' are expected
806  exposureBBox : `lsst.geom.Box2I`
807  Bounding box on which to select rows that overlap
808  exposureWcs : `lsst.afw.geom.SkyWcs`
809  World coordinate system to convert sky coords in ref cat to
810  pixel coords with which to compare with exposureBBox
811 
812  Returns
813  -------
814  refCat : `lsst.afw.table.SourceTable`
815  Source Catalog with minimal schema that overlaps exposureBBox
816  """
817  df = pd.concat(dfList)
818  # translate ra/decl coords in dataframe to detector pixel coords
819  # to down select rows that overlap the detector bbox
820  mapping = exposureWcs.getTransform().getMapping()
821  x, y = mapping.applyInverse(np.array(df[['ra', 'decl']].values*2*np.pi/360).T)
822  inBBox = lsst.geom.Box2D(exposureBBox).contains(x, y)
823  refCat = self.df2SourceCat(df[inBBox])
824  return refCat
825 
826  def df2SourceCat(self, df):
827  """Create minimal schema SourceCatalog from a pandas DataFrame.
828 
829  The forced measurement subtask expects this as input.
830 
831  Parameters
832  ----------
833  df : `pandas.DataFrame`
834  DiaObjects with locations and ids.
835 
836  Returns
837  -------
838  outputCatalog : `lsst.afw.table.SourceTable`
839  Output catalog with minimal schema.
840  """
842  outputCatalog = lsst.afw.table.SourceCatalog(schema)
843  outputCatalog.reserve(len(df))
844 
845  for diaObjectId, ra, decl in df[['ra', 'decl']].itertuples():
846  outputRecord = outputCatalog.addNew()
847  outputRecord.setId(diaObjectId)
848  outputRecord.setCoord(lsst.geom.SpherePoint(ra, decl, lsst.geom.degrees))
849  return outputCatalog
static Key< RecordId > getParentKey()
Key for the parent ID.
Definition: Source.h:273
static Schema makeMinimalSchema()
Return a minimal schema for Source tables and records.
Definition: Source.h:258
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
ConvexPolygon is a closed convex polygon on the unit sphere.
Definition: ConvexPolygon.h:57
static ConvexPolygon convexHull(std::vector< UnitVector3d > const &points)
convexHull returns the convex hull of the given set of points if it exists and throws an exception ot...
Definition: ConvexPolygon.h:65
daf::base::PropertySet * set
Definition: fits.cc:912
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:521
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
lsst::geom::Box2I bboxFromMetadata(daf::base::PropertySet &metadata)
Determine the image bounding box from its metadata (FITS header)
Definition: Image.cc:682
def run(self, coaddExposures, bbox, wcs)
Definition: getTemplate.py:603
def imageOverlapsTract(tract, imageWcs, imageBox)
def attachFootprints(self, sources, refCat, exposure, refWcs, dataRef)
def writeOutput(self, dataRef, sources)
def fetchReferences(self, dataRef, exposure)