Coverage for python / lsst / pipe / tasks / postprocess.py: 25%
816 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-07 01:55 -0800
1# This file is part of pipe_tasks.
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/>.
22__all__ = ["WriteObjectTableConfig", "WriteObjectTableTask",
23 "WriteSourceTableConfig", "WriteSourceTableTask",
24 "WriteRecalibratedSourceTableConfig", "WriteRecalibratedSourceTableTask",
25 "PostprocessAnalysis",
26 "TransformCatalogBaseConfig", "TransformCatalogBaseTask",
27 "TransformObjectCatalogConfig", "TransformObjectCatalogTask",
28 "ConsolidateObjectTableConfig", "ConsolidateObjectTableTask",
29 "TransformSourceTableConfig", "TransformSourceTableTask",
30 "ConsolidateVisitSummaryConfig", "ConsolidateVisitSummaryTask",
31 "ConsolidateSourceTableConfig", "ConsolidateSourceTableTask",
32 "MakeCcdVisitTableConfig", "MakeCcdVisitTableTask",
33 "MakeVisitTableConfig", "MakeVisitTableTask",
34 "WriteForcedSourceTableConfig", "WriteForcedSourceTableTask",
35 "TransformForcedSourceTableConfig", "TransformForcedSourceTableTask",
36 "ConsolidateTractConfig", "ConsolidateTractTask",
37 "ComputeColumnsAction", "ModelExtendednessColumnAction",
38 ]
40from collections import defaultdict
41import dataclasses
42from deprecated.sphinx import deprecated
43import functools
44import logging
45import numbers
46import os
47from typing import Iterable
49import numpy as np
50import pandas as pd
51import astropy.table
52from numpy.typing import NDArray
54import lsst.geom
55import lsst.pex.config as pexConfig
56import lsst.pipe.base as pipeBase
57import lsst.daf.base as dafBase
58from lsst.daf.butler.formatters.parquet import pandas_to_astropy
59from lsst.pex.config.configurableActions import ConfigurableAction, ConfigurableActionStructField
60from lsst.pipe.base import NoWorkFound, UpstreamFailureNoWorkFound, connectionTypes
61import lsst.afw.table as afwTable
62from lsst.afw.image import ExposureSummaryStats, ExposureF
63from lsst.meas.base import SingleFrameMeasurementTask, DetectorVisitIdGeneratorConfig
64from lsst.obs.base.utils import strip_provenance_from_fits_header, TableVStack
66from .coaddBase import reorderRefs
67from .functors import CompositeFunctor, Column
68from .schemaUtils import convertDataFrameToSdmSchema, readSdmSchemaFile
70log = logging.getLogger(__name__)
73def flattenFilters(df, noDupCols=["coord_ra", "coord_dec"], camelCase=False, inputBands=None):
74 """Flattens a dataframe with multilevel column index.
75 """
76 newDf = pd.DataFrame()
77 # band is the level 0 index
78 dfBands = df.columns.unique(level=0).values
79 for band in dfBands:
80 subdf = df[band]
81 columnFormat = "{0}{1}" if camelCase else "{0}_{1}"
82 newColumns = {c: columnFormat.format(band, c)
83 for c in subdf.columns if c not in noDupCols}
84 cols = list(newColumns.keys())
85 newDf = pd.concat([newDf, subdf[cols].rename(columns=newColumns)], axis=1)
87 # Band must be present in the input and output or else column is all NaN:
88 presentBands = dfBands if inputBands is None else list(set(inputBands).intersection(dfBands))
89 # Get the unexploded columns from any present band's partition
90 noDupDf = df[presentBands[0]][noDupCols]
91 newDf = pd.concat([noDupDf, newDf], axis=1)
92 return newDf
95class WriteObjectTableConnections(pipeBase.PipelineTaskConnections,
96 defaultTemplates={"coaddName": "deep"},
97 dimensions=("tract", "patch", "skymap")):
98 inputCatalogMeas = connectionTypes.Input(
99 doc="Catalog of source measurements on the deepCoadd.",
100 dimensions=("tract", "patch", "band", "skymap"),
101 storageClass="SourceCatalog",
102 name="{coaddName}Coadd_meas",
103 multiple=True
104 )
105 inputCatalogForcedSrc = connectionTypes.Input(
106 doc="Catalog of forced measurements (shape and position parameters held fixed) on the deepCoadd.",
107 dimensions=("tract", "patch", "band", "skymap"),
108 storageClass="SourceCatalog",
109 name="{coaddName}Coadd_forced_src",
110 multiple=True
111 )
112 inputCatalogPsfsMultiprofit = connectionTypes.Input(
113 doc="Catalog of Gaussian mixture model fit parameters for the PSF model at each object centroid.",
114 dimensions=("tract", "patch", "band", "skymap"),
115 storageClass="ArrowAstropy",
116 name="{coaddName}Coadd_psfs_multiprofit",
117 multiple=True,
118 )
119 outputCatalog = connectionTypes.Output(
120 doc="A vertical concatenation of the deepCoadd_{ref|meas|forced_src} catalogs, "
121 "stored as a DataFrame with a multi-level column index per-patch.",
122 dimensions=("tract", "patch", "skymap"),
123 storageClass="DataFrame",
124 name="{coaddName}Coadd_obj"
125 )
128class WriteObjectTableConfig(pipeBase.PipelineTaskConfig,
129 pipelineConnections=WriteObjectTableConnections):
130 coaddName = pexConfig.Field(
131 dtype=str,
132 default="deep",
133 doc="Name of coadd"
134 )
137class WriteObjectTableTask(pipeBase.PipelineTask):
138 """Write filter-merged object tables as a DataFrame in parquet format.
139 """
140 _DefaultName = "writeObjectTable"
141 ConfigClass = WriteObjectTableConfig
143 # Tag of output dataset written by `MergeSourcesTask.write`
144 outputDataset = "obj"
146 def runQuantum(self, butlerQC, inputRefs, outputRefs):
147 inputs = butlerQC.get(inputRefs)
149 catalogs = defaultdict(dict)
150 for dataset, connection in (
151 ("meas", "inputCatalogMeas"),
152 ("forced_src", "inputCatalogForcedSrc"),
153 ("psfs_multiprofit", "inputCatalogPsfsMultiprofit"),
154 ):
155 for ref, cat in zip(getattr(inputRefs, connection), inputs[connection]):
156 catalogs[ref.dataId["band"]][dataset] = cat
158 dataId = butlerQC.quantum.dataId
159 df = self.run(catalogs=catalogs, tract=dataId["tract"], patch=dataId["patch"])
160 outputs = pipeBase.Struct(outputCatalog=df)
161 butlerQC.put(outputs, outputRefs)
163 def run(self, catalogs, tract, patch):
164 """Merge multiple catalogs.
166 Parameters
167 ----------
168 catalogs : `dict`
169 Mapping from filter names to dict of catalogs.
170 tract : int
171 tractId to use for the tractId column.
172 patch : str
173 patchId to use for the patchId column.
175 Returns
176 -------
177 catalog : `pandas.DataFrame`
178 Merged dataframe.
180 Raises
181 ------
182 ValueError
183 Raised if any of the catalogs is of an unsupported type.
184 """
185 dfs = []
186 for filt, tableDict in catalogs.items():
187 for dataset, table in tableDict.items():
188 # Convert afwTable to pandas DataFrame if needed
189 if isinstance(table, pd.DataFrame):
190 df = table
191 elif isinstance(table, afwTable.SourceCatalog):
192 df = table.asAstropy().to_pandas()
193 elif isinstance(table, astropy.table.Table):
194 df = table.to_pandas()
195 else:
196 raise ValueError(f"{dataset=} has unsupported {type(table)=}")
197 df.set_index("id", drop=True, inplace=True)
199 # Sort columns by name, to ensure matching schema among patches
200 df = df.reindex(sorted(df.columns), axis=1)
201 df = df.assign(tractId=tract, patchId=patch)
203 # Make columns a 3-level MultiIndex
204 df.columns = pd.MultiIndex.from_tuples([(dataset, filt, c) for c in df.columns],
205 names=("dataset", "band", "column"))
206 dfs.append(df)
208 # We do this dance and not `pd.concat(dfs)` because the pandas
209 # concatenation uses infinite memory.
210 catalog = functools.reduce(lambda d1, d2: d1.join(d2), dfs)
211 return catalog
214class WriteSourceTableConnections(pipeBase.PipelineTaskConnections,
215 defaultTemplates={"catalogType": ""},
216 dimensions=("instrument", "visit", "detector")):
218 catalog = connectionTypes.Input(
219 doc="Input full-depth catalog of sources produced by CalibrateTask",
220 name="{catalogType}src",
221 storageClass="SourceCatalog",
222 dimensions=("instrument", "visit", "detector")
223 )
224 outputCatalog = connectionTypes.Output(
225 doc="Catalog of sources, `src` in Astropy/Parquet format. Columns are unchanged.",
226 name="{catalogType}source",
227 storageClass="ArrowAstropy",
228 dimensions=("instrument", "visit", "detector")
229 )
232class WriteSourceTableConfig(pipeBase.PipelineTaskConfig,
233 pipelineConnections=WriteSourceTableConnections):
234 pass
237class WriteSourceTableTask(pipeBase.PipelineTask):
238 """Write source table to DataFrame Parquet format.
239 """
240 _DefaultName = "writeSourceTable"
241 ConfigClass = WriteSourceTableConfig
243 def runQuantum(self, butlerQC, inputRefs, outputRefs):
244 inputs = butlerQC.get(inputRefs)
245 inputs["visit"] = butlerQC.quantum.dataId["visit"]
246 inputs["detector"] = butlerQC.quantum.dataId["detector"]
247 result = self.run(**inputs)
248 outputs = pipeBase.Struct(outputCatalog=result.table)
249 butlerQC.put(outputs, outputRefs)
251 def run(self, catalog, visit, detector, **kwargs):
252 """Convert `src` catalog to an Astropy table.
254 Parameters
255 ----------
256 catalog: `afwTable.SourceCatalog`
257 catalog to be converted
258 visit, detector: `int`
259 Visit and detector ids to be added as columns.
260 **kwargs
261 Additional keyword arguments are ignored as a convenience for
262 subclasses that pass the same arguments to several different
263 methods.
265 Returns
266 -------
267 result : `~lsst.pipe.base.Struct`
268 ``table``
269 `astropy.table.Table` version of the input catalog
270 """
271 self.log.info("Generating DataFrame from src catalog visit,detector=%i,%i", visit, detector)
272 tbl = catalog.asAstropy()
273 tbl["visit"] = visit
274 # int16 instead of uint8 because databases don't like unsigned bytes.
275 tbl["detector"] = np.int16(detector)
277 return pipeBase.Struct(table=tbl)
280class WriteRecalibratedSourceTableConnections(WriteSourceTableConnections,
281 defaultTemplates={"catalogType": ""},
282 dimensions=("instrument", "visit", "detector", "skymap")):
283 visitSummary = connectionTypes.Input(
284 doc="Input visit-summary catalog with updated calibration objects.",
285 name="finalVisitSummary",
286 storageClass="ExposureCatalog",
287 dimensions=("instrument", "visit",),
288 )
290 def __init__(self, config):
291 # We don't want the input catalog here to be an initial existence
292 # constraint in QG generation, because that can unfortunately limit the
293 # set of data IDs of inputs to other tasks, even those that run earlier
294 # (e.g. updateVisitSummary), when the input 'src' catalog is not
295 # produced. It's safer to just use 'visitSummary' existence as an
296 # initial constraint, and then let the graph prune out the detectors
297 # that don't have a 'src' for this task only.
298 self.catalog = dataclasses.replace(self.catalog, deferGraphConstraint=True)
301class WriteRecalibratedSourceTableConfig(WriteSourceTableConfig,
302 pipelineConnections=WriteRecalibratedSourceTableConnections):
304 doReevaluatePhotoCalib = pexConfig.Field(
305 dtype=bool,
306 default=True,
307 doc=("Add or replace local photoCalib columns"),
308 )
309 doReevaluateSkyWcs = pexConfig.Field(
310 dtype=bool,
311 default=True,
312 doc=("Add or replace local WCS columns and update the coord columns, coord_ra and coord_dec"),
313 )
316class WriteRecalibratedSourceTableTask(WriteSourceTableTask):
317 """Write source table to DataFrame Parquet format.
318 """
319 _DefaultName = "writeRecalibratedSourceTable"
320 ConfigClass = WriteRecalibratedSourceTableConfig
322 def runQuantum(self, butlerQC, inputRefs, outputRefs):
323 inputs = butlerQC.get(inputRefs)
325 inputs["visit"] = butlerQC.quantum.dataId["visit"]
326 inputs["detector"] = butlerQC.quantum.dataId["detector"]
328 if self.config.doReevaluatePhotoCalib or self.config.doReevaluateSkyWcs:
329 exposure = ExposureF()
330 inputs["exposure"] = self.prepareCalibratedExposure(
331 exposure=exposure,
332 visitSummary=inputs["visitSummary"],
333 detectorId=butlerQC.quantum.dataId["detector"]
334 )
335 inputs["catalog"] = self.addCalibColumns(**inputs)
337 result = self.run(**inputs)
338 outputs = pipeBase.Struct(outputCatalog=result.table)
339 butlerQC.put(outputs, outputRefs)
341 def prepareCalibratedExposure(self, exposure, detectorId, visitSummary=None):
342 """Prepare a calibrated exposure and apply external calibrations
343 if so configured.
345 Parameters
346 ----------
347 exposure : `lsst.afw.image.exposure.Exposure`
348 Input exposure to adjust calibrations. May be an empty Exposure.
349 detectorId : `int`
350 Detector ID associated with the exposure.
351 visitSummary : `lsst.afw.table.ExposureCatalog`, optional
352 Exposure catalog with all calibration objects. WCS and PhotoCalib
353 are always applied if ``visitSummary`` is provided and those
354 components are not `None`.
356 Returns
357 -------
358 exposure : `lsst.afw.image.exposure.Exposure`
359 Exposure with adjusted calibrations.
360 """
361 if visitSummary is not None:
362 row = visitSummary.find(detectorId)
363 if row is None:
364 raise pipeBase.NoWorkFound(f"Visit summary for detector {detectorId} is missing.")
365 if (photoCalib := row.getPhotoCalib()) is None:
366 self.log.warning("Detector id %s has None for photoCalib in visit summary; "
367 "skipping reevaluation of photoCalib.", detectorId)
368 exposure.setPhotoCalib(None)
369 else:
370 exposure.setPhotoCalib(photoCalib)
371 if (skyWcs := row.getWcs()) is None:
372 self.log.warning("Detector id %s has None for skyWcs in visit summary; "
373 "skipping reevaluation of skyWcs.", detectorId)
374 exposure.setWcs(None)
375 else:
376 exposure.setWcs(skyWcs)
378 return exposure
380 def addCalibColumns(self, catalog, exposure, **kwargs):
381 """Add replace columns with calibs evaluated at each centroid
383 Add or replace 'base_LocalWcs' and 'base_LocalPhotoCalib' columns in
384 a source catalog, by rerunning the plugins.
386 Parameters
387 ----------
388 catalog : `lsst.afw.table.SourceCatalog`
389 catalog to which calib columns will be added
390 exposure : `lsst.afw.image.exposure.Exposure`
391 Exposure with attached PhotoCalibs and SkyWcs attributes to be
392 reevaluated at local centroids. Pixels are not required.
393 **kwargs
394 Additional keyword arguments are ignored to facilitate passing the
395 same arguments to several methods.
397 Returns
398 -------
399 newCat: `lsst.afw.table.SourceCatalog`
400 Source Catalog with requested local calib columns
401 """
402 measureConfig = SingleFrameMeasurementTask.ConfigClass()
403 measureConfig.doReplaceWithNoise = False
405 # Clear all slots, because we aren't running the relevant plugins.
406 for slot in measureConfig.slots:
407 setattr(measureConfig.slots, slot, None)
409 measureConfig.plugins.names = []
410 if self.config.doReevaluateSkyWcs:
411 measureConfig.plugins.names.add("base_LocalWcs")
412 self.log.info("Re-evaluating base_LocalWcs plugin")
413 if self.config.doReevaluatePhotoCalib:
414 measureConfig.plugins.names.add("base_LocalPhotoCalib")
415 self.log.info("Re-evaluating base_LocalPhotoCalib plugin")
416 pluginsNotToCopy = tuple(measureConfig.plugins.names)
418 # Create a new schema and catalog
419 # Copy all columns from original except for the ones to reevaluate
420 aliasMap = catalog.schema.getAliasMap()
421 mapper = afwTable.SchemaMapper(catalog.schema)
422 for item in catalog.schema:
423 if not item.field.getName().startswith(pluginsNotToCopy):
424 mapper.addMapping(item.key)
426 schema = mapper.getOutputSchema()
427 measurement = SingleFrameMeasurementTask(config=measureConfig, schema=schema)
428 schema.setAliasMap(aliasMap)
429 newCat = afwTable.SourceCatalog(schema)
430 newCat.extend(catalog, mapper=mapper)
432 # Fluxes in sourceCatalogs are in counts, so there are no fluxes to
433 # update here. LocalPhotoCalibs are applied during transform tasks.
434 # Update coord_ra/coord_dec, which are expected to be positions on the
435 # sky and are used as such in sdm tables without transform
436 if self.config.doReevaluateSkyWcs and exposure.wcs is not None:
437 afwTable.updateSourceCoords(exposure.wcs, newCat)
438 wcsPlugin = measurement.plugins["base_LocalWcs"]
439 else:
440 wcsPlugin = None
442 if self.config.doReevaluatePhotoCalib and exposure.getPhotoCalib() is not None:
443 pcPlugin = measurement.plugins["base_LocalPhotoCalib"]
444 else:
445 pcPlugin = None
447 for row in newCat:
448 if wcsPlugin is not None:
449 wcsPlugin.measure(row, exposure)
450 if pcPlugin is not None:
451 pcPlugin.measure(row, exposure)
453 return newCat
456class PostprocessAnalysis(object):
457 """Calculate columns from DataFrames or handles storing DataFrames.
459 This object manages and organizes an arbitrary set of computations
460 on a catalog. The catalog is defined by a
461 `DeferredDatasetHandle` or `InMemoryDatasetHandle` object
462 (or list thereof), such as a ``deepCoadd_obj`` dataset, and the
463 computations are defined by a collection of
464 `~lsst.pipe.tasks.functors.Functor` objects (or, equivalently, a
465 ``CompositeFunctor``).
467 After the object is initialized, accessing the ``.df`` attribute (which
468 holds the `pandas.DataFrame` containing the results of the calculations)
469 triggers computation of said dataframe.
471 One of the conveniences of using this object is the ability to define a
472 desired common filter for all functors. This enables the same functor
473 collection to be passed to several different `PostprocessAnalysis` objects
474 without having to change the original functor collection, since the ``filt``
475 keyword argument of this object triggers an overwrite of the ``filt``
476 property for all functors in the collection.
478 This object also allows a list of refFlags to be passed, and defines a set
479 of default refFlags that are always included even if not requested.
481 If a list of DataFrames or Handles is passed, rather than a single one,
482 then the calculations will be mapped over all the input catalogs. In
483 principle, it should be straightforward to parallelize this activity, but
484 initial tests have failed (see TODO in code comments).
486 Parameters
487 ----------
488 handles : `~lsst.daf.butler.DeferredDatasetHandle` or
489 `~lsst.pipe.base.InMemoryDatasetHandle` or
490 list of these.
491 Source catalog(s) for computation.
492 functors : `list`, `dict`, or `~lsst.pipe.tasks.functors.CompositeFunctor`
493 Computations to do (functors that act on ``handles``).
494 If a dict, the output
495 DataFrame will have columns keyed accordingly.
496 If a list, the column keys will come from the
497 ``.shortname`` attribute of each functor.
499 filt : `str`, optional
500 Filter in which to calculate. If provided,
501 this will overwrite any existing ``.filt`` attribute
502 of the provided functors.
504 flags : `list`, optional
505 List of flags (per-band) to include in output table.
506 Taken from the ``meas`` dataset if applied to a multilevel Object Table.
508 refFlags : `list`, optional
509 List of refFlags (only reference band) to include in output table.
511 forcedFlags : `list`, optional
512 List of flags (per-band) to include in output table.
513 Taken from the ``forced_src`` dataset if applied to a
514 multilevel Object Table. Intended for flags from measurement plugins
515 only run during multi-band forced-photometry.
516 """
517 _defaultRefFlags = []
518 _defaultFuncs = ()
520 def __init__(self, handles, functors, filt=None, flags=None, refFlags=None, forcedFlags=None):
521 self.handles = handles
522 self.functors = functors
524 self.filt = filt
525 self.flags = list(flags) if flags is not None else []
526 self.forcedFlags = list(forcedFlags) if forcedFlags is not None else []
527 self.refFlags = list(self._defaultRefFlags)
528 if refFlags is not None:
529 self.refFlags += list(refFlags)
531 self._df = None
533 @property
534 def defaultFuncs(self):
535 funcs = dict(self._defaultFuncs)
536 return funcs
538 @property
539 def func(self):
540 additionalFuncs = self.defaultFuncs
541 additionalFuncs.update({flag: Column(flag, dataset="forced_src") for flag in self.forcedFlags})
542 additionalFuncs.update({flag: Column(flag, dataset="ref") for flag in self.refFlags})
543 additionalFuncs.update({flag: Column(flag, dataset="meas") for flag in self.flags})
545 if isinstance(self.functors, CompositeFunctor):
546 func = self.functors
547 else:
548 func = CompositeFunctor(self.functors)
550 func.funcDict.update(additionalFuncs)
551 func.filt = self.filt
553 return func
555 @property
556 def noDupCols(self):
557 return [name for name, func in self.func.funcDict.items() if func.noDup]
559 @property
560 def df(self):
561 if self._df is None:
562 self.compute()
563 return self._df
565 def compute(self, dropna=False, pool=None):
566 # map over multiple handles
567 if type(self.handles) in (list, tuple):
568 if pool is None:
569 dflist = [self.func(handle, dropna=dropna) for handle in self.handles]
570 else:
571 # TODO: Figure out why this doesn't work (pyarrow pickling
572 # issues?)
573 dflist = pool.map(functools.partial(self.func, dropna=dropna), self.handles)
574 self._df = pd.concat(dflist)
575 else:
576 self._df = self.func(self.handles, dropna=dropna)
578 return self._df
581class TransformCatalogBaseConnections(pipeBase.PipelineTaskConnections,
582 dimensions=()):
583 """Expected Connections for subclasses of TransformCatalogBaseTask.
585 Must be subclassed.
586 """
587 inputCatalog = connectionTypes.Input(
588 name="",
589 storageClass="DataFrame",
590 )
591 outputCatalog = connectionTypes.Output(
592 name="",
593 storageClass="ArrowAstropy",
594 )
597class TransformCatalogBaseConfig(pipeBase.PipelineTaskConfig,
598 pipelineConnections=TransformCatalogBaseConnections):
599 functorFile = pexConfig.Field(
600 dtype=str,
601 doc="Path to YAML file specifying Science Data Model functors to use "
602 "when copying columns and computing calibrated values.",
603 default=None,
604 optional=True
605 )
606 primaryKey = pexConfig.Field(
607 dtype=str,
608 doc="Name of column to be set as the DataFrame index. If None, the index"
609 "will be named `id`",
610 default=None,
611 optional=True
612 )
613 columnsFromDataId = pexConfig.ListField(
614 dtype=str,
615 default=None,
616 optional=True,
617 doc="Columns to extract from the dataId",
618 )
621class TransformCatalogBaseTask(pipeBase.PipelineTask):
622 """Base class for transforming/standardizing a catalog by applying functors
623 that convert units and apply calibrations.
625 The purpose of this task is to perform a set of computations on an input
626 ``DeferredDatasetHandle`` or ``InMemoryDatasetHandle`` that holds a
627 ``DataFrame`` dataset (such as ``deepCoadd_obj``), and write the results to
628 a new dataset (which needs to be declared in an ``outputDataset``
629 attribute).
631 The calculations to be performed are defined in a YAML file that specifies
632 a set of functors to be computed, provided as a ``--functorFile`` config
633 parameter. An example of such a YAML file is the following:
635 funcs:
636 sourceId:
637 functor: Index
638 x:
639 functor: Column
640 args: slot_Centroid_x
641 y:
642 functor: Column
643 args: slot_Centroid_y
644 psfFlux:
645 functor: LocalNanojansky
646 args:
647 - slot_PsfFlux_instFlux
648 - slot_PsfFlux_instFluxErr
649 - base_LocalPhotoCalib
650 - base_LocalPhotoCalibErr
651 psfFluxErr:
652 functor: LocalNanojanskyErr
653 args:
654 - slot_PsfFlux_instFlux
655 - slot_PsfFlux_instFluxErr
656 - base_LocalPhotoCalib
657 - base_LocalPhotoCalibErr
658 flags:
659 - detect_isPrimary
661 The names for each entry under "func" will become the names of columns in
662 the output dataset. All the functors referenced are defined in
663 `~lsst.pipe.tasks.functors`. Positional arguments to be passed to each
664 functor are in the `args` list, and any additional entries for each column
665 other than "functor" or "args" (e.g., ``'filt'``, ``'dataset'``) are
666 treated as keyword arguments to be passed to the functor initialization.
668 The "flags" entry is the default shortcut for `Column` functors.
669 All columns listed under "flags" will be copied to the output table
670 untransformed. They can be of any datatype.
671 In the special case of transforming a multi-level oject table with
672 band and dataset indices (deepCoadd_obj), these will be taked from the
673 ``meas`` dataset and exploded out per band.
675 There are two special shortcuts that only apply when transforming
676 multi-level Object (deepCoadd_obj) tables:
677 - The "refFlags" entry is shortcut for `Column` functor
678 taken from the ``ref`` dataset if transforming an ObjectTable.
679 - The "forcedFlags" entry is shortcut for `Column` functors.
680 taken from the ``forced_src`` dataset if transforming an ObjectTable.
681 These are expanded out per band.
684 This task uses the `lsst.pipe.tasks.postprocess.PostprocessAnalysis` object
685 to organize and excecute the calculations.
686 """
687 @property
688 def _DefaultName(self):
689 raise NotImplementedError("Subclass must define the \"_DefaultName\" attribute.")
691 @property
692 def outputDataset(self):
693 raise NotImplementedError("Subclass must define the \"outputDataset\" attribute.")
695 @property
696 def inputDataset(self):
697 raise NotImplementedError("Subclass must define \"inputDataset\" attribute.")
699 @property
700 def ConfigClass(self):
701 raise NotImplementedError("Subclass must define \"ConfigClass\" attribute.")
703 def __init__(self, *args, **kwargs):
704 super().__init__(*args, **kwargs)
705 if self.config.functorFile:
706 self.log.info("Loading tranform functor definitions from %s",
707 self.config.functorFile)
708 self.funcs = CompositeFunctor.from_file(self.config.functorFile)
709 self.funcs.update(dict(PostprocessAnalysis._defaultFuncs))
710 else:
711 self.funcs = None
713 def runQuantum(self, butlerQC, inputRefs, outputRefs):
714 inputs = butlerQC.get(inputRefs)
715 if self.funcs is None:
716 raise ValueError("config.functorFile is None. "
717 "Must be a valid path to yaml in order to run Task as a PipelineTask.")
718 result = self.run(handle=inputs["inputCatalog"], funcs=self.funcs,
719 dataId=dict(outputRefs.outputCatalog.dataId.mapping))
720 butlerQC.put(result, outputRefs)
722 def run(self, handle, funcs=None, dataId=None, band=None):
723 """Do postprocessing calculations
725 Takes a ``DeferredDatasetHandle`` or ``InMemoryDatasetHandle`` or
726 ``DataFrame`` object and dataId,
727 returns a dataframe with results of postprocessing calculations.
729 Parameters
730 ----------
731 handles : `~lsst.daf.butler.DeferredDatasetHandle` or
732 `~lsst.pipe.base.InMemoryDatasetHandle` or
733 `~pandas.DataFrame`, or list of these.
734 DataFrames from which calculations are done.
735 funcs : `~lsst.pipe.tasks.functors.Functor`
736 Functors to apply to the table's columns
737 dataId : dict, optional
738 Used to add a `patchId` column to the output dataframe.
739 band : `str`, optional
740 Filter band that is being processed.
742 Returns
743 -------
744 result : `lsst.pipe.base.Struct`
745 Result struct, with a single ``outputCatalog`` attribute holding
746 the transformed catalog.
747 """
748 self.log.info("Transforming/standardizing the source table dataId: %s", dataId)
750 df = self.transform(band, handle, funcs, dataId).df
751 self.log.info("Made a table of %d columns and %d rows", len(df.columns), len(df))
753 if len(df) == 0:
754 raise UpstreamFailureNoWorkFound(
755 "Input catalog is empty, so there is nothing to transform/standardize",
756 )
758 result = pipeBase.Struct(outputCatalog=pandas_to_astropy(df))
759 return result
761 def getFunctors(self):
762 return self.funcs
764 def getAnalysis(self, handles, funcs=None, band=None):
765 if funcs is None:
766 funcs = self.funcs
767 analysis = PostprocessAnalysis(handles, funcs, filt=band)
768 return analysis
770 def transform(self, band, handles, funcs, dataId):
771 analysis = self.getAnalysis(handles, funcs=funcs, band=band)
772 df = analysis.df
773 if dataId and self.config.columnsFromDataId:
774 for key in self.config.columnsFromDataId:
775 if key in dataId:
776 if key == "detector":
777 # int16 instead of uint8 because databases don't like unsigned bytes.
778 df[key] = np.int16(dataId[key])
779 else:
780 df[key] = dataId[key]
781 else:
782 raise ValueError(f"'{key}' in config.columnsFromDataId not found in dataId: {dataId}")
784 if self.config.primaryKey:
785 if df.index.name != self.config.primaryKey and self.config.primaryKey in df:
786 df.reset_index(inplace=True, drop=True)
787 df.set_index(self.config.primaryKey, inplace=True)
789 return pipeBase.Struct(
790 df=df,
791 analysis=analysis
792 )
795class TransformObjectCatalogConnections(pipeBase.PipelineTaskConnections,
796 defaultTemplates={"coaddName": "deep"},
797 dimensions=("tract", "patch", "skymap")):
798 inputCatalog = connectionTypes.Input(
799 doc="The vertical concatenation of the {coaddName}_{meas|forced_src|psfs_multiprofit} catalogs, "
800 "stored as a DataFrame with a multi-level column index per-patch.",
801 dimensions=("tract", "patch", "skymap"),
802 storageClass="DataFrame",
803 name="{coaddName}Coadd_obj",
804 deferLoad=True,
805 )
806 inputCatalogRef = connectionTypes.Input(
807 doc="Catalog marking the primary detection (which band provides a good shape and position)"
808 "for each detection in deepCoadd_mergeDet.",
809 dimensions=("tract", "patch", "skymap"),
810 storageClass="SourceCatalog",
811 name="{coaddName}Coadd_ref",
812 deferLoad=True,
813 )
814 inputCatalogExpMultiprofit = connectionTypes.Input(
815 doc="Catalog of multiband Exponential fits.",
816 dimensions=("tract", "patch", "skymap"),
817 storageClass="ArrowAstropy",
818 name="{coaddName}Coadd_Exp_multiprofit",
819 deferLoad=True,
820 )
821 inputCatalogSersicMultiprofit = connectionTypes.Input(
822 doc="Catalog of multiband Sersic fits.",
823 dimensions=("tract", "patch", "skymap"),
824 storageClass="ArrowAstropy",
825 name="{coaddName}Coadd_Sersic_multiprofit",
826 deferLoad=True,
827 )
828 inputCatalogEpoch = connectionTypes.Input(
829 doc="Catalog of mean epochs for each object per band.",
830 dimensions=("tract", "patch", "skymap"),
831 storageClass="ArrowAstropy",
832 name="object_epoch",
833 deferLoad=True,
834 )
835 outputCatalog = connectionTypes.Output(
836 doc="Per-Patch Object Table of columns transformed from the deepCoadd_obj table per the standard "
837 "data model.",
838 dimensions=("tract", "patch", "skymap"),
839 storageClass="ArrowAstropy",
840 name="objectTable"
841 )
843 def __init__(self, *, config=None):
844 super().__init__(config=config)
845 if config.multilevelOutput:
846 self.outputCatalog = dataclasses.replace(self.outputCatalog, storageClass="DataFrame")
849class TransformObjectCatalogConfig(TransformCatalogBaseConfig,
850 pipelineConnections=TransformObjectCatalogConnections):
851 coaddName = pexConfig.Field(
852 dtype=str,
853 default="deep",
854 doc="Name of coadd"
855 )
856 outputBands = pexConfig.ListField(
857 dtype=str,
858 default=None,
859 optional=True,
860 doc=("These bands and only these bands will appear in the output,"
861 " NaN-filled if the input does not include them."
862 " If None, then use all bands found in the input.")
863 )
864 camelCase = pexConfig.Field(
865 dtype=bool,
866 default=False,
867 doc=("Write per-band columns names with camelCase, else underscore "
868 "For example: gPsFlux instead of g_PsFlux.")
869 )
870 multilevelOutput = pexConfig.Field(
871 dtype=bool,
872 default=False,
873 doc=("Whether results dataframe should have a multilevel column index (True) or be flat "
874 "and name-munged (False). If True, the output storage class will be "
875 "set to DataFrame, since astropy tables do not support multi-level indexing."),
876 deprecated="Support for multi-level outputs is deprecated and will be removed after v29.",
877 )
878 goodFlags = pexConfig.ListField(
879 dtype=str,
880 default=[],
881 doc=("List of 'good' flags that should be set False when populating empty tables. "
882 "All other flags are considered to be 'bad' flags and will be set to True.")
883 )
884 floatFillValue = pexConfig.Field(
885 dtype=float,
886 default=np.nan,
887 doc="Fill value for float fields when populating empty tables."
888 )
889 integerFillValue = pexConfig.Field(
890 dtype=int,
891 default=-1,
892 doc="Fill value for integer fields when populating empty tables."
893 )
895 def setDefaults(self):
896 super().setDefaults()
897 self.functorFile = os.path.join("$PIPE_TASKS_DIR", "schemas", "Object.yaml")
898 self.primaryKey = "objectId"
899 self.columnsFromDataId = ["tract", "patch"]
900 self.goodFlags = ["calib_astrometry_used",
901 "calib_photometry_reserved",
902 "calib_photometry_used",
903 "calib_psf_candidate",
904 "calib_psf_reserved",
905 "calib_psf_used"]
908class TransformObjectCatalogTask(TransformCatalogBaseTask):
909 """Produce a flattened Object Table to match the format specified in
910 sdm_schemas.
912 Do the same set of postprocessing calculations on all bands.
914 This is identical to `TransformCatalogBaseTask`, except for that it does
915 the specified functor calculations for all filters present in the
916 input `deepCoadd_obj` table. Any specific ``"filt"`` keywords specified
917 by the YAML file will be superceded.
918 """
919 _DefaultName = "transformObjectCatalog"
920 ConfigClass = TransformObjectCatalogConfig
922 # ref must go first because other datasets may need columns from it
923 datasets_multiband = ("ref", "epoch", "Exp_multiprofit", "Sersic_multiprofit")
925 def runQuantum(self, butlerQC, inputRefs, outputRefs):
926 inputs = butlerQC.get(inputRefs)
927 if self.funcs is None:
928 raise ValueError("config.functorFile is None. "
929 "Must be a valid path to yaml in order to run Task as a PipelineTask.")
930 result = self.run(handle=inputs["inputCatalog"], funcs=self.funcs,
931 dataId=dict(outputRefs.outputCatalog.dataId.mapping),
932 handle_ref=inputs["inputCatalogRef"],
933 handle_epoch=inputs["inputCatalogEpoch"],
934 handle_Exp_multiprofit=inputs["inputCatalogExpMultiprofit"],
935 handle_Sersic_multiprofit=inputs["inputCatalogSersicMultiprofit"],
936 )
937 butlerQC.put(result, outputRefs)
939 def run(self, handle, funcs=None, dataId=None, band=None, **kwargs):
940 # NOTE: band kwarg is ignored here.
941 # TODO: Document and improve funcs argument usage in DM-48895
942 # self.getAnalysis only supports list, dict and CompositeFunctor
943 if isinstance(funcs, CompositeFunctor):
944 funcDict_in = funcs.funcDict
945 elif isinstance(funcs, dict):
946 funcDict_in = funcs
947 elif isinstance(funcs, list):
948 funcDict_in = {idx: v for idx, v in enumerate(funcs)}
949 else:
950 raise TypeError(f"Unsupported {type(funcs)=}")
952 handles_multi = {}
953 funcDicts_multiband = {}
954 for dataset in self.datasets_multiband:
955 if (handle_multi := kwargs.get(f"handle_{dataset}")) is None:
956 raise RuntimeError(f"Missing required handle_{dataset} kwarg")
957 handles_multi[dataset] = handle_multi
958 funcDicts_multiband[dataset] = {}
960 dfDict = {}
961 analysisDict = {}
962 templateDf = pd.DataFrame()
964 columns = handle.get(component="columns")
965 inputBands = columns.unique(level=1).values
967 outputBands = self.config.outputBands if self.config.outputBands else inputBands
969 # Split up funcs for per-band and multiband tables
970 funcDict_band = {}
971 # Check if dataset needs to add any meas table columns
972 columns_add_from_ref = defaultdict(list)
974 for name, func in funcDict_in.items():
975 if func.dataset in funcDicts_multiband:
976 # This is something like a MultibandColumn
977 if band := getattr(func, "band_to_check", None):
978 if band not in outputBands:
979 continue
980 # This is something like a ReferenceBand that has configurable bands
981 elif hasattr(func, "bands"):
982 # TODO: Determine if this can be avoided DM-48895
983 # This will work fine if the init doesn't manipulate bands
984 # If it does, then one would need to make a new functor
985 # Determining the (kw)args is tricky in that case
986 func.bands = tuple(inputBands)
987 if hasattr(func, "columns_ref"):
988 columns_add_from_ref[func.dataset].extend(func.columns_ref)
990 funcDict = funcDicts_multiband.get(func.dataset, funcDict_band)
991 funcDict[name] = func
993 funcs_band = CompositeFunctor(funcDict_band)
995 # Perform transform for data of filters that exist in the handle dataframe.
996 for inputBand in inputBands:
997 if inputBand not in outputBands:
998 self.log.info("Ignoring %s band data in the input", inputBand)
999 continue
1000 self.log.info("Transforming the catalog of band %s", inputBand)
1001 result = self.transform(inputBand, handle, funcs_band, dataId)
1002 dfDict[inputBand] = result.df
1003 analysisDict[inputBand] = result.analysis
1004 if templateDf.empty:
1005 templateDf = result.df
1007 # Put filler values in columns of other wanted bands
1008 for filt in outputBands:
1009 if filt not in dfDict:
1010 self.log.info("Adding empty columns for band %s", filt)
1011 dfTemp = templateDf.copy()
1012 for col in dfTemp.columns:
1013 testValue = dfTemp[col].values[0]
1014 if isinstance(testValue, (np.bool_, pd.BooleanDtype)):
1015 # Boolean flag type, check if it is a "good" flag
1016 if col in self.config.goodFlags:
1017 fillValue = False
1018 else:
1019 fillValue = True
1020 elif isinstance(testValue, numbers.Integral):
1021 # Checking numbers.Integral catches all flavors
1022 # of python, numpy, pandas, etc. integers.
1023 # We must ensure this is not an unsigned integer.
1024 if isinstance(testValue, np.unsignedinteger):
1025 raise ValueError("Parquet tables may not have unsigned integer columns.")
1026 else:
1027 fillValue = self.config.integerFillValue
1028 else:
1029 fillValue = self.config.floatFillValue
1030 dfTemp[col].values[:] = fillValue
1031 dfDict[filt] = dfTemp
1033 # This makes a multilevel column index, with band as first level
1034 df = pd.concat(dfDict, axis=1, names=["band", "column"])
1035 name_index = df.index.name
1037 # TODO: Remove in DM-48895
1038 if not self.config.multilevelOutput:
1039 noDupCols = list(set.union(*[set(v.noDupCols) for v in analysisDict.values()]))
1040 if self.config.primaryKey in noDupCols:
1041 noDupCols.remove(self.config.primaryKey)
1042 if dataId and self.config.columnsFromDataId:
1043 noDupCols += self.config.columnsFromDataId
1044 df = flattenFilters(df, noDupCols=noDupCols, camelCase=self.config.camelCase,
1045 inputBands=inputBands)
1047 # This will be assigned by running ref first
1048 df_ref = None
1050 # Apply per-dataset functors to each multiband dataset in turn
1051 for dataset, funcDict in funcDicts_multiband.items():
1052 handle_multiband = handles_multi[dataset]
1053 df_dataset = handle_multiband.get()
1054 if isinstance(df_dataset, astropy.table.Table):
1055 # Allow astropy table inputs to already have the output index
1056 if name_index not in df_dataset.colnames:
1057 if self.config.primaryKey in df_dataset.colnames:
1058 name_index_ap = self.config.primaryKey
1059 else:
1060 raise RuntimeError(
1061 f"Neither of {name_index=} nor {self.config.primaryKey=} appear in"
1062 f" {df_dataset.colnames=} for {dataset=}"
1063 )
1064 else:
1065 name_index_ap = name_index
1066 df_dataset = df_dataset.to_pandas().set_index(name_index_ap, drop=False)
1067 elif isinstance(df_dataset, afwTable.SourceCatalog):
1068 df_dataset = df_dataset.asAstropy().to_pandas().set_index(name_index, drop=False)
1070 if dataset == "ref":
1071 df_ref = df_dataset
1072 else:
1073 if df_ref is None:
1074 raise RuntimeError(f"ref must be the first dataset, not {dataset}")
1075 # Add columns as if they're a sorted set
1076 for column in {key: None for key in columns_add_from_ref.get(dataset, [])}:
1077 df_dataset[column] = df_ref[column]
1079 # TODO: should funcDict have noDup funcs removed?
1080 # noDup was intended for per-band tables.
1081 result = self.transform(
1082 None,
1083 pipeBase.InMemoryDatasetHandle(df_dataset, storageClass="DataFrame"),
1084 CompositeFunctor(funcDict),
1085 dataId,
1086 )
1087 result.df.index.name = name_index
1088 # Drop columns from dataId if present (patch, tract)
1089 if self.config.columnsFromDataId:
1090 columns_drop = [column for column in self.config.columnsFromDataId if column in result.df]
1091 if columns_drop:
1092 result.df.drop(columns_drop, axis=1, inplace=True)
1093 # Make the same multi-index for the multiband table if needed
1094 # This might end up making copies, one of several reasons to avoid
1095 # using multilevel indexes, or DataFrames at all
1096 to_concat = pd.concat(
1097 {band: result.df for band in self.config.outputBands}, axis=1, names=["band", "column"]
1098 ) if self.config.multilevelOutput else result.df
1099 df = pd.concat([df, to_concat], axis=1)
1100 analysisDict[dataset] = result.analysis
1101 del result
1103 df.index.name = self.config.primaryKey
1105 if not self.config.multilevelOutput:
1106 tbl = pandas_to_astropy(df)
1107 else:
1108 tbl = df
1110 self.log.info("Made a table of %d columns and %d rows", len(tbl.columns), len(tbl))
1112 return pipeBase.Struct(outputCatalog=tbl)
1115class ComputeColumnsAction(ConfigurableAction):
1116 """An action that computes multiple vectors from an input.
1118 This class is meant to be compatible with analysis_tools'
1119 AnalysisAction class, which cannot be a dependency of pipe_tasks."""
1121 def getInputSchema(self) -> dict[str, type[NDArray]]:
1122 """Return the required inputs for this action.
1124 This function is meant to be compatible with
1125 """
1126 raise NotImplementedError("This method must be overloaded in subclasses")
1128 def __call__(self, table: astropy.table.Table) -> dict[str, NDArray]:
1129 """This method must return a dict of computed columns."""
1130 raise NotImplementedError("This method must be overloaded in subclasses")
1133class ExtendednessColumnActionBase(ComputeColumnsAction):
1134 bands = pexConfig.ListField[str](
1135 doc="The bands to make single-band outputs for.",
1136 default=["u", "g", "r", "i", "z", "y"]
1137 )
1138 bands_combined = pexConfig.DictField[str, str](
1139 doc="Multiband classification column specialization. Keys specify the"
1140 " name of the column and values are a comma-separated list of"
1141 " bands, all of which must be contained in the bands listed.",
1142 default={"griz": "g,r,i,z"},
1143 itemCheck=lambda x: (len(y := x.split(",")) > 1) & (len(set(y)) == len(y)),
1144 )
1145 model_column_flux = pexConfig.Field[str](
1146 doc="The model flux column to use for computing the difference to"
1147 " to the S/N flux. Must contain the {band} and {model} templates.",
1148 default="{band}_{model}Flux",
1149 check=lambda x: ("{band}" in x) and ("{model}" in x,),
1150 )
1151 model_column_flux_err = pexConfig.Field[str](
1152 doc="The model flux error column to use for computing the difference"
1153 " to the S/N flux. Must contain the {band} and {model} templates.",
1154 default="{band}_{model}FluxErr",
1155 check=lambda x: ("{band}" in x) and ("{model}" in x,),
1156 )
1157 model_flux_name = pexConfig.Field[str](
1158 doc="The extended object model to use to compared to PSF model fluxes",
1159 default="sersic",
1160 )
1161 output_column = pexConfig.Field[str](
1162 doc="Name of the output column. Must contain the {band} template",
1163 default="{band}_model_extendedness",
1164 check=lambda x: "{band}" in x,
1165 )
1166 psf_column_flux = pexConfig.Field[str](
1167 doc="The name of the PSF flux column. Must contain the {band} template.",
1168 default="{band}_psfFlux",
1169 check=lambda x: "{band}" in x,
1170 )
1171 psf_column_flux_err = pexConfig.Field[str](
1172 doc="The name of the PSF flux error column. Must contain the {band} template.",
1173 default="{band}_psfFluxErr",
1174 check=lambda x: "{band}" in x,
1175 )
1176 size_column = pexConfig.Field[str](
1177 doc="The column to use for applying size cuts. Must contain the {axis} template.",
1178 default="exponential_reff_{axis}",
1179 )
1181 def getInputSchema(self) -> Iterable[tuple[str, type[NDArray]]]:
1182 size_column = self.size_column
1183 schema = [
1184 (size_column.format(axis=axis), NDArray[float]) for axis in ("x", "y")
1185 ]
1186 model = self.model_flux_name
1187 for column in (
1188 self.psf_column_flux, self.psf_column_flux_err,
1189 self.model_column_flux, self.model_column_flux_err,
1190 ):
1191 schema.extend([
1192 (column.format(band=band, model=model), NDArray[float]) for band in self.bands
1193 ])
1195 return schema
1197 def validate(self):
1198 super().validate()
1199 errors = []
1200 for name, band_combined in self.bands_combined.items():
1201 bands = band_combined.split(",")
1202 bands_missing = [band for band in bands if band not in self.bands]
1203 if bands_missing:
1204 errors.append(
1205 f"self.bands_combined[{name}] contains bands={bands_missing} not in {self.bands=}"
1206 )
1207 if errors:
1208 raise ValueError(f"Validation failed due to errors: {'; '.join(errors)}")
1210 def _get_fluxes(self, table, band: str):
1211 model = self.model_flux_name
1212 flux_psf, fluxerr_psf, flux_model, fluxerr_model = (
1213 np.array(table[column.format(band=band, model=model)])
1214 for column in (
1215 self.psf_column_flux, self.psf_column_flux_err,
1216 self.model_column_flux, self.model_column_flux_err,
1217 )
1218 )
1219 return flux_psf, fluxerr_psf, flux_model, fluxerr_model
1222class ModelExtendednessColumnAction(ExtendednessColumnActionBase):
1223 fluxerr_coefficent = pexConfig.Field[float](
1224 doc="The coefficient to multiply the flux error by when adding to the model flux.",
1225 default=0.5,
1226 check=lambda x: x >= 0,
1227 )
1228 fluxerr_stretch = pexConfig.Field[float](
1229 doc="The factor to multiply flux error-scaled ratios by to derive extendedness.",
1230 default=5.0,
1231 check=lambda x: x > 0,
1232 )
1233 good_sn_min = pexConfig.Field[float](
1234 doc="Minimum PSF S/N to include objects if"
1235 " min_n_good_to_shift_flux_ratio is > 0; ignored otherwise.",
1236 default=10.,
1237 )
1238 max_reff_compact = pexConfig.Field[float](
1239 doc="The maximum effective radius in pixels below which an object is"
1240 " classified as not extended, regardless of other parameter values.",
1241 default=0.25,
1242 )
1243 min_n_good_to_shift_flux_ratio = pexConfig.Field[int](
1244 doc="Minimum number of objects with PSF S/N > good_sn_min and with "
1245 " size larger than max_reff_compact to use to compute the median "
1246 " PSF-to-model flux ratio, which is assumed to be 1 otherwise."
1247 " If this value is not >0, the median flux ratio will be kept 1.",
1248 default=0,
1249 )
1251 def __call__(self, table: astropy.table.Table) -> dict[str, NDArray]:
1252 size_column = self.size_column
1253 size_model = np.sqrt(
1254 0.5*(table[size_column.format(axis='x')]**2 + table[size_column.format(axis='y')]**2)
1255 )
1256 small = size_model < self.max_reff_compact
1257 n_obj = len(table)
1258 band_mappings_to_process = {band: [band] for band in self.bands}
1259 band_mappings_to_process.update({k: v.split(",") for k, v in self.bands_combined.items()})
1260 output = {}
1261 for output_band, input_bands in band_mappings_to_process.items():
1262 if len(input_bands) > 1:
1263 flux_psf, fluxerr_psf_sq, flux_model, fluxerr_model_sq = (
1264 np.zeros(n_obj, dtype=float) for _ in range(4))
1265 for input_band in input_bands:
1266 flux_psf_b, fluxerr_psf_b, flux_model_b, fluxerr_model_b = self._get_fluxes(
1267 table, band=input_band)
1268 # There's no point adding S/N < 0 fluxes
1269 good = np.isfinite(flux_psf_b) & np.isfinite(flux_model_b) & (
1270 flux_psf_b > 0) & (flux_model_b > 0) & (fluxerr_psf_b > 0) & (fluxerr_model_b > 0)
1271 flux_psf[good] += flux_psf_b[good]
1272 flux_model[good] += flux_model_b[good]
1273 fluxerr_psf_sq[good] += fluxerr_psf_b[good]**2
1274 fluxerr_model_sq[good] += fluxerr_model_b[good]**2
1275 fluxerr_psf = np.sqrt(fluxerr_psf_sq)
1276 fluxerr_model = np.sqrt(fluxerr_model_sq)
1277 fluxerr_psf[fluxerr_psf == 0] = np.inf
1278 fluxerr_model[fluxerr_model == 0] = np.inf
1279 else:
1280 flux_psf, fluxerr_psf, flux_model, fluxerr_model = self._get_fluxes(
1281 table, band=input_bands[0])
1283 psf_sn = flux_psf/fluxerr_psf
1284 flux_ratio = np.array(flux_psf / flux_model)
1286 if self.min_n_good_to_shift_flux_ratio > 0:
1287 good = small & (psf_sn > self.good_sn_min)
1288 # Attempt to correct any flux-independent systematic offset
1289 # Might need to be a function of S/N
1290 if np.sum(good == True) > self.min_n_good_to_shift_flux_ratio: # noqa: E712
1291 flux_ratio *= 1./np.nanmedian(flux_ratio[good])
1293 flux_ratio_err = np.sqrt(
1294 (fluxerr_psf/flux_model)**2 + (fluxerr_model*fluxerr_psf/flux_model**2)**2
1295 )
1296 extendedness = (1 - flux_ratio) + self.fluxerr_coefficent*flux_ratio_err
1297 extendedness *= np.sqrt(size_model/self.max_reff_compact)
1298 extendedness[(extendedness < 0) & (extendedness > -np.inf)] = 0
1299 # Make it sigmoid-like with a stretch
1300 stretch = self.fluxerr_stretch
1301 extendedness *= stretch
1302 extendedness = np.clip((stretch + 1)/stretch*extendedness/(1 + extendedness), 0, 1)
1304 column_out = self.output_column.format(band=output_band)
1305 output[column_out] = extendedness
1306 return output
1309class ConsolidateObjectTableConnections(pipeBase.PipelineTaskConnections,
1310 dimensions=("tract", "skymap")):
1311 inputCatalogs = connectionTypes.Input(
1312 doc="Per-Patch objectTables conforming to the standard data model.",
1313 name="objectTable",
1314 storageClass="ArrowAstropy",
1315 dimensions=("tract", "patch", "skymap"),
1316 multiple=True,
1317 deferLoad=True,
1318 )
1319 outputCatalog = connectionTypes.Output(
1320 doc="Pre-tract horizontal concatenation of the input objectTables",
1321 name="objectTable_tract",
1322 storageClass="ArrowAstropy",
1323 dimensions=("tract", "skymap"),
1324 )
1327class ConsolidateObjectTableConfig(pipeBase.PipelineTaskConfig,
1328 pipelineConnections=ConsolidateObjectTableConnections):
1329 actions = ConfigurableActionStructField[ComputeColumnsAction](
1330 doc="Actions to add columns to the final object table",
1331 )
1332 coaddName = pexConfig.Field[str](
1333 default="deep",
1334 doc="Name of coadd"
1335 )
1337 def setDefaults(self):
1338 super().setDefaults()
1339 self.actions.extendedness = ModelExtendednessColumnAction()
1342class ConsolidateObjectTableTask(pipeBase.PipelineTask):
1343 """Write patch-merged source tables to a tract-level DataFrame Parquet file.
1345 Concatenates `objectTable` list into a per-visit `objectTable_tract`.
1346 """
1347 _DefaultName = "consolidateObjectTable"
1348 ConfigClass = ConsolidateObjectTableConfig
1350 inputDataset = "objectTable"
1351 outputDataset = "objectTable_tract"
1353 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1354 inputs = butlerQC.get(inputRefs)
1355 self.log.info("Concatenating %s per-patch Object Tables",
1356 len(inputs["inputCatalogs"]))
1357 table = TableVStack.vstack_handles(inputs["inputCatalogs"])
1358 for action in self.config.actions:
1359 computed = action(table)
1360 for key, values in computed.items():
1361 table[key] = values.astype(np.float32) if values.dtype == np.float64 else values
1362 butlerQC.put(pipeBase.Struct(outputCatalog=table), outputRefs)
1365class TransformSourceTableConnections(pipeBase.PipelineTaskConnections,
1366 defaultTemplates={"catalogType": ""},
1367 dimensions=("instrument", "visit", "detector")):
1369 inputCatalog = connectionTypes.Input(
1370 doc="Wide input catalog of sources produced by WriteSourceTableTask",
1371 name="{catalogType}source",
1372 storageClass="DataFrame",
1373 dimensions=("instrument", "visit", "detector"),
1374 deferLoad=True
1375 )
1376 outputCatalog = connectionTypes.Output(
1377 doc="Narrower, per-detector Source Table transformed and converted per a "
1378 "specified set of functors",
1379 name="{catalogType}sourceTable",
1380 storageClass="ArrowAstropy",
1381 dimensions=("instrument", "visit", "detector")
1382 )
1385class TransformSourceTableConfig(TransformCatalogBaseConfig,
1386 pipelineConnections=TransformSourceTableConnections):
1388 def setDefaults(self):
1389 super().setDefaults()
1390 self.functorFile = os.path.join("$PIPE_TASKS_DIR", "schemas", "Source.yaml")
1391 self.primaryKey = "sourceId"
1392 self.columnsFromDataId = ["visit", "detector", "band", "physical_filter"]
1395class TransformSourceTableTask(TransformCatalogBaseTask):
1396 """Transform/standardize a source catalog
1397 """
1398 _DefaultName = "transformSourceTable"
1399 ConfigClass = TransformSourceTableConfig
1402class ConsolidateVisitSummaryConnections(pipeBase.PipelineTaskConnections,
1403 dimensions=("instrument", "visit",),
1404 defaultTemplates={"calexpType": ""}):
1405 calexp = connectionTypes.Input(
1406 doc="Processed exposures used for metadata",
1407 name="calexp",
1408 storageClass="ExposureF",
1409 dimensions=("instrument", "visit", "detector"),
1410 deferLoad=True,
1411 multiple=True,
1412 )
1413 visitSummary = connectionTypes.Output(
1414 doc=("Per-visit consolidated exposure metadata. These catalogs use "
1415 "detector id for the id and are sorted for fast lookups of a "
1416 "detector."),
1417 name="visitSummary",
1418 storageClass="ExposureCatalog",
1419 dimensions=("instrument", "visit"),
1420 )
1421 visitSummarySchema = connectionTypes.InitOutput(
1422 doc="Schema of the visitSummary catalog",
1423 name="visitSummary_schema",
1424 storageClass="ExposureCatalog",
1425 )
1428class ConsolidateVisitSummaryConfig(pipeBase.PipelineTaskConfig,
1429 pipelineConnections=ConsolidateVisitSummaryConnections):
1430 """Config for ConsolidateVisitSummaryTask"""
1432 full = pexConfig.Field(
1433 "Whether to propate all exposure components. "
1434 "This adds PSF, aperture correction map, transmission curve, and detector, which can increase file "
1435 "size by more than factor of 10, but it makes the visit summaries produced by this task fully usable"
1436 "by tasks that were designed to run downstream of lsst.drp.tasks.UpdateVisitSummaryTask.",
1437 dtype=bool,
1438 default=False,
1439 )
1442class ConsolidateVisitSummaryTask(pipeBase.PipelineTask):
1443 """Task to consolidate per-detector visit metadata.
1445 This task aggregates the following metadata from all the detectors in a
1446 single visit into an exposure catalog:
1447 - The visitInfo.
1448 - The wcs.
1449 - The photoCalib.
1450 - The physical_filter and band (if available).
1451 - The PSF model.
1452 - The aperture correction map.
1453 - The transmission curve.
1454 - The psf size, shape, and effective area at the center of the detector.
1455 - The corners of the bounding box in right ascension/declination.
1457 Tests for this task are performed in ci_hsc_gen3.
1458 """
1459 _DefaultName = "consolidateVisitSummary"
1460 ConfigClass = ConsolidateVisitSummaryConfig
1462 def __init__(self, **kwargs):
1463 super().__init__(**kwargs)
1464 self.schema = afwTable.ExposureTable.makeMinimalSchema()
1465 self.schema.addField("visit", type="L", doc="Visit number")
1466 self.schema.addField("physical_filter", type="String", size=32, doc="Physical filter")
1467 self.schema.addField("band", type="String", size=32, doc="Name of band")
1468 ExposureSummaryStats.update_schema(self.schema)
1469 self.visitSummarySchema = afwTable.ExposureCatalog(self.schema)
1471 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1472 dataRefs = butlerQC.get(inputRefs.calexp)
1473 visit = dataRefs[0].dataId["visit"]
1475 self.log.debug("Concatenating metadata from %d per-detector calexps (visit %d)",
1476 len(dataRefs), visit)
1478 expCatalog = self._combineExposureMetadata(visit, dataRefs)
1480 butlerQC.put(expCatalog, outputRefs.visitSummary)
1482 def _combineExposureMetadata(self, visit, dataRefs):
1483 """Make a combined exposure catalog from a list of dataRefs.
1484 These dataRefs must point to exposures with wcs, summaryStats,
1485 and other visit metadata.
1487 Parameters
1488 ----------
1489 visit : `int`
1490 Visit identification number.
1491 dataRefs : `list` of `lsst.daf.butler.DeferredDatasetHandle`
1492 List of dataRefs in visit.
1494 Returns
1495 -------
1496 visitSummary : `lsst.afw.table.ExposureCatalog`
1497 Exposure catalog with per-detector summary information.
1498 """
1499 cat = afwTable.ExposureCatalog(self.schema)
1500 cat.resize(len(dataRefs))
1502 cat["visit"] = visit
1504 for i, dataRef in enumerate(dataRefs):
1505 visitInfo = dataRef.get(component="visitInfo")
1506 filterLabel = dataRef.get(component="filter")
1507 summaryStats = dataRef.get(component="summaryStats")
1508 detector = dataRef.get(component="detector")
1509 wcs = dataRef.get(component="wcs")
1510 photoCalib = dataRef.get(component="photoCalib")
1511 bbox = dataRef.get(component="bbox")
1512 validPolygon = dataRef.get(component="validPolygon")
1514 rec = cat[i]
1515 rec.setBBox(bbox)
1516 rec.setVisitInfo(visitInfo)
1517 rec.setWcs(wcs)
1518 rec.setPhotoCalib(photoCalib)
1519 rec.setValidPolygon(validPolygon)
1521 if self.config.full:
1522 rec.setPsf(dataRef.get(component="psf"))
1523 rec.setApCorrMap(dataRef.get(component="apCorrMap"))
1524 rec.setTransmissionCurve(dataRef.get(component="transmissionCurve"))
1526 rec["physical_filter"] = filterLabel.physicalLabel if filterLabel.hasPhysicalLabel() else ""
1527 rec["band"] = filterLabel.bandLabel if filterLabel.hasBandLabel() else ""
1528 rec.setId(detector.getId())
1529 summaryStats.update_record(rec)
1531 if not cat:
1532 raise pipeBase.NoWorkFound(
1533 "No detectors had sufficient information to make a visit summary row."
1534 )
1536 metadata = dafBase.PropertyList()
1537 metadata.add("COMMENT", "Catalog id is detector id, sorted.")
1538 # We are looping over existing datarefs, so the following is true
1539 metadata.add("COMMENT", "Only detectors with data have entries.")
1540 cat.setMetadata(metadata)
1542 cat.sort()
1543 return cat
1546class ConsolidateSourceTableConnections(pipeBase.PipelineTaskConnections,
1547 defaultTemplates={"catalogType": ""},
1548 dimensions=("instrument", "visit")):
1549 inputCatalogs = connectionTypes.Input(
1550 doc="Input per-detector Source Tables",
1551 name="{catalogType}sourceTable",
1552 storageClass="ArrowAstropy",
1553 dimensions=("instrument", "visit", "detector"),
1554 multiple=True,
1555 deferLoad=True,
1556 )
1557 outputCatalog = connectionTypes.Output(
1558 doc="Per-visit concatenation of Source Table",
1559 name="{catalogType}sourceTable_visit",
1560 storageClass="ArrowAstropy",
1561 dimensions=("instrument", "visit")
1562 )
1565class ConsolidateSourceTableConfig(pipeBase.PipelineTaskConfig,
1566 pipelineConnections=ConsolidateSourceTableConnections):
1567 pass
1570class ConsolidateSourceTableTask(pipeBase.PipelineTask):
1571 """Concatenate `sourceTable` list into a per-visit `sourceTable_visit`
1572 """
1573 _DefaultName = "consolidateSourceTable"
1574 ConfigClass = ConsolidateSourceTableConfig
1576 inputDataset = "sourceTable"
1577 outputDataset = "sourceTable_visit"
1579 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1580 # Docstring inherited.
1581 detectorOrder = [ref.dataId["detector"] for ref in inputRefs.inputCatalogs]
1582 detectorOrder.sort()
1583 inputRefs = reorderRefs(inputRefs, detectorOrder, dataIdKey="detector")
1584 inputs = butlerQC.get(inputRefs)
1585 self.log.info("Concatenating %s per-detector Source Tables",
1586 len(inputs["inputCatalogs"]))
1587 table = TableVStack.vstack_handles(inputs["inputCatalogs"])
1588 butlerQC.put(pipeBase.Struct(outputCatalog=table), outputRefs)
1591class MakeCcdVisitTableConnections(pipeBase.PipelineTaskConnections,
1592 dimensions=("instrument",),
1593 defaultTemplates={"calexpType": ""}):
1594 visitSummaryRefs = connectionTypes.Input(
1595 doc="Data references for per-visit consolidated exposure metadata",
1596 name="finalVisitSummary",
1597 storageClass="ExposureCatalog",
1598 dimensions=("instrument", "visit"),
1599 multiple=True,
1600 deferLoad=True,
1601 )
1602 outputCatalog = connectionTypes.Output(
1603 doc="CCD and Visit metadata table",
1604 name="ccdVisitTable",
1605 storageClass="ArrowAstropy",
1606 dimensions=("instrument",)
1607 )
1610class MakeCcdVisitTableConfig(pipeBase.PipelineTaskConfig,
1611 pipelineConnections=MakeCcdVisitTableConnections):
1612 idGenerator = DetectorVisitIdGeneratorConfig.make_field()
1615class MakeCcdVisitTableTask(pipeBase.PipelineTask):
1616 """Produce a `ccdVisitTable` from the visit summary exposure catalogs.
1617 """
1618 _DefaultName = "makeCcdVisitTable"
1619 ConfigClass = MakeCcdVisitTableConfig
1621 def run(self, visitSummaryRefs):
1622 """Make a table of ccd information from the visit summary catalogs.
1624 Parameters
1625 ----------
1626 visitSummaryRefs : `list` of `lsst.daf.butler.DeferredDatasetHandle`
1627 List of DeferredDatasetHandles pointing to exposure catalogs with
1628 per-detector summary information.
1630 Returns
1631 -------
1632 result : `~lsst.pipe.base.Struct`
1633 Results struct with attribute:
1635 ``outputCatalog``
1636 Catalog of ccd and visit information.
1637 """
1638 ccdEntries = []
1639 for visitSummaryRef in visitSummaryRefs:
1640 visitSummary = visitSummaryRef.get()
1641 if not visitSummary:
1642 continue
1643 visitInfo = visitSummary[0].getVisitInfo()
1645 # Strip provenance to prevent merge confusion.
1646 strip_provenance_from_fits_header(visitSummary.metadata)
1648 ccdEntry = {}
1649 summaryTable = visitSummary.asAstropy()
1650 selectColumns = ["id", "visit", "physical_filter", "band", "ra", "dec",
1651 "pixelScale", "zenithDistance",
1652 "expTime", "zeroPoint", "psfSigma", "skyBg", "skyNoise",
1653 "astromOffsetMean", "astromOffsetStd", "nPsfStar",
1654 "psfStarDeltaE1Median", "psfStarDeltaE2Median",
1655 "psfStarDeltaE1Scatter", "psfStarDeltaE2Scatter",
1656 "psfStarDeltaSizeMedian", "psfStarDeltaSizeScatter",
1657 "psfStarScaledDeltaSizeScatter", "psfTraceRadiusDelta",
1658 "psfApFluxDelta", "psfApCorrSigmaScaledDelta",
1659 "maxDistToNearestPsf", "starEMedian", "starUnNormalizedEMedian",
1660 "effTime", "effTimePsfSigmaScale",
1661 "effTimeSkyBgScale", "effTimeZeroPointScale",
1662 "magLim"]
1663 ccdEntry = summaryTable[selectColumns]
1664 # 'visit' is the human readable visit number.
1665 # 'visitId' is the key to the visitId table. They are the same.
1666 # Technically you should join to get the visit from the visit
1667 # table.
1668 ccdEntry.rename_column("visit", "visitId")
1669 ccdEntry.rename_column("id", "detectorId")
1671 # RFC-924: Temporarily keep a duplicate "decl" entry for backwards
1672 # compatibility. To be removed after September 2023.
1673 ccdEntry["decl"] = ccdEntry["dec"]
1675 ccdEntry["ccdVisitId"] = [
1676 self.config.idGenerator.apply(
1677 visitSummaryRef.dataId,
1678 detector=detector_id,
1679 is_exposure=False,
1680 ).catalog_id # The "catalog ID" here is the ccdVisit ID
1681 # because it's usually the ID for a whole catalog
1682 # with a {visit, detector}, and that's the main
1683 # use case for IdGenerator. This usage for a
1684 # summary table is rare.
1685 for detector_id in summaryTable["id"]
1686 ]
1687 ccdEntry["detector"] = summaryTable["id"]
1688 ccdEntry["seeing"] = (
1689 visitSummary["psfSigma"] * visitSummary["pixelScale"] * np.sqrt(8 * np.log(2))
1690 )
1691 ccdEntry["skyRotation"] = visitInfo.getBoresightRotAngle().asDegrees()
1692 ccdEntry["expMidpt"] = np.datetime64(visitInfo.date.nsecs(scale=dafBase.DateTime.TAI), "ns")
1693 ccdEntry["expMidptMJD"] = visitInfo.getDate().get(dafBase.DateTime.MJD)
1694 expTime = visitInfo.getExposureTime()
1695 ccdEntry["obsStart"] = (
1696 ccdEntry["expMidpt"] - 0.5 * np.timedelta64(int(expTime * 1E9), "ns")
1697 )
1698 expTime_days = expTime / (60*60*24)
1699 ccdEntry["obsStartMJD"] = ccdEntry["expMidptMJD"] - 0.5 * expTime_days
1700 ccdEntry["darkTime"] = visitInfo.getDarkTime()
1701 ccdEntry["xSize"] = summaryTable["bbox_max_x"] - summaryTable["bbox_min_x"]
1702 ccdEntry["ySize"] = summaryTable["bbox_max_y"] - summaryTable["bbox_min_y"]
1703 ccdEntry["llcra"] = summaryTable["raCorners"][:, 0]
1704 ccdEntry["llcdec"] = summaryTable["decCorners"][:, 0]
1705 ccdEntry["ulcra"] = summaryTable["raCorners"][:, 1]
1706 ccdEntry["ulcdec"] = summaryTable["decCorners"][:, 1]
1707 ccdEntry["urcra"] = summaryTable["raCorners"][:, 2]
1708 ccdEntry["urcdec"] = summaryTable["decCorners"][:, 2]
1709 ccdEntry["lrcra"] = summaryTable["raCorners"][:, 3]
1710 ccdEntry["lrcdec"] = summaryTable["decCorners"][:, 3]
1711 # TODO: DM-30618, Add raftName, nExposures, ccdTemp, binX, binY,
1712 # and flags, and decide if WCS, and llcx, llcy, ulcx, ulcy, etc.
1713 # values are actually wanted.
1714 ccdEntries.append(ccdEntry)
1716 outputCatalog = astropy.table.vstack(ccdEntries, join_type="exact")
1717 return pipeBase.Struct(outputCatalog=outputCatalog)
1720class MakeVisitTableConnections(pipeBase.PipelineTaskConnections,
1721 dimensions=("instrument",),
1722 defaultTemplates={"calexpType": ""}):
1723 visitSummaries = connectionTypes.Input(
1724 doc="Per-visit consolidated exposure metadata",
1725 name="finalVisitSummary",
1726 storageClass="ExposureCatalog",
1727 dimensions=("instrument", "visit",),
1728 multiple=True,
1729 deferLoad=True,
1730 )
1731 outputCatalog = connectionTypes.Output(
1732 doc="Visit metadata table",
1733 name="visitTable",
1734 storageClass="ArrowAstropy",
1735 dimensions=("instrument",)
1736 )
1739class MakeVisitTableConfig(pipeBase.PipelineTaskConfig,
1740 pipelineConnections=MakeVisitTableConnections):
1741 pass
1744class MakeVisitTableTask(pipeBase.PipelineTask):
1745 """Produce a `visitTable` from the visit summary exposure catalogs.
1746 """
1747 _DefaultName = "makeVisitTable"
1748 ConfigClass = MakeVisitTableConfig
1750 def run(self, visitSummaries):
1751 """Make a table of visit information from the visit summary catalogs.
1753 Parameters
1754 ----------
1755 visitSummaries : `list` of `lsst.afw.table.ExposureCatalog`
1756 List of exposure catalogs with per-detector summary information.
1757 Returns
1758 -------
1759 result : `~lsst.pipe.base.Struct`
1760 Results struct with attribute:
1762 ``outputCatalog``
1763 Catalog of visit information.
1764 """
1765 visitEntries = []
1766 for visitSummary in visitSummaries:
1767 visitSummary = visitSummary.get()
1768 if not visitSummary:
1769 continue
1770 visitRow = visitSummary[0]
1771 visitInfo = visitRow.getVisitInfo()
1773 visitEntry = {}
1774 visitEntry["visitId"] = visitRow["visit"]
1775 visitEntry["visit"] = visitRow["visit"]
1776 visitEntry["physical_filter"] = visitRow["physical_filter"]
1777 visitEntry["band"] = visitRow["band"]
1778 raDec = visitInfo.getBoresightRaDec()
1779 visitEntry["ra"] = raDec.getRa().asDegrees()
1780 visitEntry["dec"] = raDec.getDec().asDegrees()
1782 # RFC-924: Temporarily keep a duplicate "decl" entry for backwards
1783 # compatibility. To be removed after September 2023.
1784 visitEntry["decl"] = visitEntry["dec"]
1786 visitEntry["skyRotation"] = visitInfo.getBoresightRotAngle().asDegrees()
1787 azAlt = visitInfo.getBoresightAzAlt()
1788 visitEntry["azimuth"] = azAlt.getLongitude().asDegrees()
1789 visitEntry["altitude"] = azAlt.getLatitude().asDegrees()
1790 visitEntry["zenithDistance"] = 90 - azAlt.getLatitude().asDegrees()
1791 visitEntry["airmass"] = visitInfo.getBoresightAirmass()
1792 expTime = visitInfo.getExposureTime()
1793 visitEntry["expTime"] = expTime
1794 visitEntry["expMidpt"] = np.datetime64(visitInfo.date.nsecs(scale=dafBase.DateTime.TAI), "ns")
1795 visitEntry["expMidptMJD"] = visitInfo.getDate().get(dafBase.DateTime.MJD)
1796 visitEntry["obsStart"] = visitEntry["expMidpt"] - 0.5 * np.timedelta64(int(expTime * 1E9), "ns")
1797 expTime_days = expTime / (60*60*24)
1798 visitEntry["obsStartMJD"] = visitEntry["expMidptMJD"] - 0.5 * expTime_days
1799 visitEntries.append(visitEntry)
1801 # TODO: DM-30623, Add programId, exposureType, cameraTemp,
1802 # mirror1Temp, mirror2Temp, mirror3Temp, domeTemp, externalTemp,
1803 # dimmSeeing, pwvGPS, pwvMW, flags, nExposures.
1805 outputCatalog = astropy.table.Table(rows=visitEntries)
1806 return pipeBase.Struct(outputCatalog=outputCatalog)
1809@deprecated(reason="This task is replaced by lsst.pipe.tasks.ForcedPhotCcdTask. "
1810 "This task will be removed after v30.",
1811 version="v29.0", category=FutureWarning)
1812class WriteForcedSourceTableConnections(pipeBase.PipelineTaskConnections,
1813 dimensions=("instrument", "visit", "detector", "skymap", "tract")):
1815 inputCatalog = connectionTypes.Input(
1816 doc="Primary per-detector, single-epoch forced-photometry catalog. "
1817 "By default, it is the output of ForcedPhotCcdTask on calexps",
1818 name="forced_src",
1819 storageClass="SourceCatalog",
1820 dimensions=("instrument", "visit", "detector", "skymap", "tract")
1821 )
1822 inputCatalogDiff = connectionTypes.Input(
1823 doc="Secondary multi-epoch, per-detector, forced photometry catalog. "
1824 "By default, it is the output of ForcedPhotCcdTask run on image differences.",
1825 name="forced_diff",
1826 storageClass="SourceCatalog",
1827 dimensions=("instrument", "visit", "detector", "skymap", "tract")
1828 )
1829 outputCatalog = connectionTypes.Output(
1830 doc="InputCatalogs horizonatally joined on `objectId` in DataFrame parquet format",
1831 name="mergedForcedSource",
1832 storageClass="DataFrame",
1833 dimensions=("instrument", "visit", "detector", "skymap", "tract")
1834 )
1837@deprecated(reason="This task is replaced by lsst.pipe.tasks.ForcedPhotCcdTask. "
1838 "This task will be removed after v30.",
1839 version="v29.0", category=FutureWarning)
1840class WriteForcedSourceTableConfig(pipeBase.PipelineTaskConfig,
1841 pipelineConnections=WriteForcedSourceTableConnections):
1842 key = lsst.pex.config.Field(
1843 doc="Column on which to join the two input tables on and make the primary key of the output",
1844 dtype=str,
1845 default="objectId",
1846 )
1849@deprecated(reason="This task is replaced by lsst.pipe.tasks.ForcedPhotCcdTask. "
1850 "This task will be removed after v30.",
1851 version="v29.0", category=FutureWarning)
1852class WriteForcedSourceTableTask(pipeBase.PipelineTask):
1853 """Merge and convert per-detector forced source catalogs to DataFrame Parquet format.
1855 Because the predecessor ForcedPhotCcdTask operates per-detector,
1856 per-tract, (i.e., it has tract in its dimensions), detectors
1857 on the tract boundary may have multiple forced source catalogs.
1859 The successor task TransformForcedSourceTable runs per-patch
1860 and temporally-aggregates overlapping mergedForcedSource catalogs from all
1861 available multiple epochs.
1862 """
1863 _DefaultName = "writeForcedSourceTable"
1864 ConfigClass = WriteForcedSourceTableConfig
1866 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1867 inputs = butlerQC.get(inputRefs)
1868 inputs["visit"] = butlerQC.quantum.dataId["visit"]
1869 inputs["detector"] = butlerQC.quantum.dataId["detector"]
1870 inputs["band"] = butlerQC.quantum.dataId["band"]
1871 outputs = self.run(**inputs)
1872 butlerQC.put(outputs, outputRefs)
1874 def run(self, inputCatalog, inputCatalogDiff, visit, detector, band=None):
1875 dfs = []
1876 for table, dataset, in zip((inputCatalog, inputCatalogDiff), ("calexp", "diff")):
1877 df = table.asAstropy().to_pandas().set_index(self.config.key, drop=False)
1878 df = df.reindex(sorted(df.columns), axis=1)
1879 df["visit"] = visit
1880 # int16 instead of uint8 because databases don't like unsigned bytes.
1881 df["detector"] = np.int16(detector)
1882 df["band"] = band if band else pd.NA
1883 df.columns = pd.MultiIndex.from_tuples([(dataset, c) for c in df.columns],
1884 names=("dataset", "column"))
1886 dfs.append(df)
1888 outputCatalog = functools.reduce(lambda d1, d2: d1.join(d2), dfs)
1889 return pipeBase.Struct(outputCatalog=outputCatalog)
1892class TransformForcedSourceTableConnections(pipeBase.PipelineTaskConnections,
1893 dimensions=("instrument", "skymap", "patch", "tract")):
1895 inputCatalogs = connectionTypes.Input(
1896 doc="DataFrames of merged ForcedSources produced by WriteForcedSourceTableTask",
1897 name="mergedForcedSource",
1898 storageClass="DataFrame",
1899 dimensions=("instrument", "visit", "detector", "skymap", "tract"),
1900 multiple=True,
1901 deferLoad=True
1902 )
1903 referenceCatalog = connectionTypes.Input(
1904 doc="Reference catalog which was used to seed the forcedPhot. Columns "
1905 "objectId, detect_isPrimary, detect_isTractInner, detect_isPatchInner "
1906 "are expected.",
1907 name="objectTable",
1908 storageClass="DataFrame",
1909 dimensions=("tract", "patch", "skymap"),
1910 deferLoad=True
1911 )
1912 outputCatalog = connectionTypes.Output(
1913 doc="Narrower, temporally-aggregated, per-patch ForcedSource Table transformed and converted per a "
1914 "specified set of functors",
1915 name="forcedSourceTable",
1916 storageClass="ArrowAstropy",
1917 dimensions=("tract", "patch", "skymap")
1918 )
1921class TransformForcedSourceTableConfig(TransformCatalogBaseConfig,
1922 pipelineConnections=TransformForcedSourceTableConnections):
1923 referenceColumns = pexConfig.ListField(
1924 dtype=str,
1925 default=["detect_isPrimary", "detect_isTractInner", "detect_isPatchInner"],
1926 optional=True,
1927 doc="Columns to pull from reference catalog",
1928 )
1929 keyRef = lsst.pex.config.Field(
1930 doc="Column on which to join the two input tables on and make the primary key of the output",
1931 dtype=str,
1932 default="objectId",
1933 )
1934 key = lsst.pex.config.Field(
1935 doc="Rename the output DataFrame index to this name",
1936 dtype=str,
1937 default="forcedSourceId",
1938 )
1940 def setDefaults(self):
1941 super().setDefaults()
1942 self.functorFile = os.path.join("$PIPE_TASKS_DIR", "schemas", "ForcedSource.yaml")
1943 self.columnsFromDataId = ["tract", "patch"]
1946class TransformForcedSourceTableTask(TransformCatalogBaseTask):
1947 """Transform/standardize a ForcedSource catalog
1949 Transforms each wide, per-detector forcedSource DataFrame per the
1950 specification file (per-camera defaults found in ForcedSource.yaml).
1951 All epochs that overlap the patch are aggregated into one per-patch
1952 narrow-DataFrame file.
1954 No de-duplication of rows is performed. Duplicate resolutions flags are
1955 pulled in from the referenceCatalog: `detect_isPrimary`,
1956 `detect_isTractInner`,`detect_isPatchInner`, so that user may de-duplicate
1957 for analysis or compare duplicates for QA.
1959 The resulting table includes multiple bands. Epochs (MJDs) and other useful
1960 per-visit rows can be retreived by joining with the CcdVisitTable on
1961 ccdVisitId.
1962 """
1963 _DefaultName = "transformForcedSourceTable"
1964 ConfigClass = TransformForcedSourceTableConfig
1966 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1967 inputs = butlerQC.get(inputRefs)
1968 if self.funcs is None:
1969 raise ValueError("config.functorFile is None. "
1970 "Must be a valid path to yaml in order to run Task as a PipelineTask.")
1971 outputs = self.run(inputs["inputCatalogs"], inputs["referenceCatalog"], funcs=self.funcs,
1972 dataId=dict(outputRefs.outputCatalog.dataId.mapping))
1974 butlerQC.put(outputs, outputRefs)
1976 def run(self, inputCatalogs, referenceCatalog, funcs=None, dataId=None, band=None):
1977 dfs = []
1978 refColumns = list(self.config.referenceColumns)
1979 refColumns.append(self.config.keyRef)
1980 ref = referenceCatalog.get(parameters={"columns": refColumns})
1981 if ref.index.name != self.config.keyRef:
1982 # If the DataFrame we loaded was originally written as some other
1983 # Parquet type, it probably doesn't have the index set. If it was
1984 # written as a DataFrame, the index should already be set and
1985 # trying to set it again would be an error, since it doens't exist
1986 # as a regular column anymore.
1987 ref.set_index(self.config.keyRef, inplace=True)
1988 self.log.info("Aggregating %s input catalogs" % (len(inputCatalogs)))
1989 for handle in inputCatalogs:
1990 result = self.transform(None, handle, funcs, dataId)
1991 # Filter for only rows that were detected on (overlap) the patch
1992 dfs.append(result.df.join(ref, how="inner"))
1994 outputCatalog = pd.concat(dfs)
1996 if outputCatalog.empty:
1997 raise NoWorkFound(f"No forced photometry rows for {dataId}.")
1999 # Now that we are done joining on config.keyRef
2000 # Change index to config.key by
2001 outputCatalog.index.rename(self.config.keyRef, inplace=True)
2002 # Add config.keyRef to the column list
2003 outputCatalog.reset_index(inplace=True)
2005 if "forcedSourceId" in outputCatalog.columns:
2006 # Set the forcedSourceId to the index. This is specified in the
2007 # ForcedSource.yaml
2008 outputCatalog.set_index("forcedSourceId", inplace=True, verify_integrity=True)
2009 # Rename it to the config.key
2010 outputCatalog.index.rename(self.config.key, inplace=True)
2012 self.log.info("Made a table of %d columns and %d rows",
2013 len(outputCatalog.columns), len(outputCatalog))
2014 return pipeBase.Struct(outputCatalog=pandas_to_astropy(outputCatalog))
2017class ConsolidateTractConnections(pipeBase.PipelineTaskConnections,
2018 defaultTemplates={"catalogType": ""},
2019 dimensions=("instrument", "tract")):
2020 inputCatalogs = connectionTypes.Input(
2021 doc="Input per-patch DataFrame Tables to be concatenated",
2022 name="{catalogType}ForcedSourceTable",
2023 storageClass="DataFrame",
2024 dimensions=("tract", "patch", "skymap"),
2025 multiple=True,
2026 )
2028 outputCatalog = connectionTypes.Output(
2029 doc="Output per-tract concatenation of DataFrame Tables",
2030 name="{catalogType}ForcedSourceTable_tract",
2031 storageClass="DataFrame",
2032 dimensions=("tract", "skymap"),
2033 )
2036class ConsolidateTractConfig(pipeBase.PipelineTaskConfig,
2037 pipelineConnections=ConsolidateTractConnections):
2039 doUseSchema = pexConfig.Field(
2040 dtype=bool,
2041 default=False,
2042 doc="Use an existing schema to coerce the data types of the output columns."
2043 )
2044 schemaDir = pexConfig.Field(
2045 dtype=str,
2046 doc="Path to the directory containing schema definitions.",
2047 default=os.path.join("${SDM_SCHEMAS_DIR}",
2048 "yml"),
2049 optional=True,
2050 )
2051 schemaFile = pexConfig.Field(
2052 dtype=str,
2053 doc="Yaml file specifying the schema of the output catalog.",
2054 optional=True,
2055 )
2056 tableName = pexConfig.Field(
2057 dtype=str,
2058 doc="Name of the table in the schema file to read.",
2059 optional=True,
2060 )
2063class ConsolidateTractTask(pipeBase.PipelineTask):
2064 """Concatenate any per-patch, dataframe list into a single
2065 per-tract DataFrame.
2066 """
2067 _DefaultName = "ConsolidateTract"
2068 ConfigClass = ConsolidateTractConfig
2070 def runQuantum(self, butlerQC, inputRefs, outputRefs):
2071 inputs = butlerQC.get(inputRefs)
2072 # Not checking at least one inputCatalog exists because that'd be an
2073 # empty QG.
2074 self.log.info("Concatenating %s per-patch %s Tables",
2075 len(inputs["inputCatalogs"]),
2076 inputRefs.inputCatalogs[0].datasetType.name)
2077 df = pd.concat(inputs["inputCatalogs"])
2078 if self.config.doUseSchema:
2079 schemaFile = os.path.join(self.config.schemaDir, self.config.schemaFile)
2080 schema = readSdmSchemaFile(schemaFile)
2081 df = convertDataFrameToSdmSchema(schema, df, tableName=self.config.tableName)
2082 butlerQC.put(pipeBase.Struct(outputCatalog=df), outputRefs)
2085class ConsolidateParentTractConnections(
2086 pipeBase.PipelineTaskConnections,
2087 dimensions=("instrument", "tract")
2088):
2089 inputCatalogs = connectionTypes.Input(
2090 doc="Parents of the deblended objects",
2091 name="object_parent_patch",
2092 storageClass="SourceCatalog",
2093 dimensions=("tract", "patch", "skymap"),
2094 multiple=True,
2095 )
2097 outputCatalog = connectionTypes.Output(
2098 doc="Output per-tract concatenation of DataFrame Tables",
2099 name="object_parent",
2100 storageClass="ArrowAstropy",
2101 dimensions=("tract", "skymap"),
2102 )
2105class ConsolidateParentTractConfig(
2106 pipeBase.PipelineTaskConfig,
2107 pipelineConnections=ConsolidateParentTractConnections,
2108):
2109 pass
2112class ConsolidateParentTractTask(pipeBase.PipelineTask):
2113 """Concatenate any per-patch, dataframe list into a single
2114 per-tract DataFrame.
2115 """
2116 _DefaultName = "ConsolidateTract"
2117 ConfigClass = ConsolidateParentTractConfig
2119 def runQuantum(self, butlerQC, inputRefs, outputRefs):
2120 self.log.info("Concatenating %s per-patch %s Tables",
2121 len(inputRefs.inputCatalogs),
2122 inputRefs.inputCatalogs[0].datasetType.name)
2124 tables = []
2125 for ref in inputRefs.inputCatalogs:
2126 catalog = butlerQC.get(ref)
2127 table = catalog.asAstropy()
2129 # Rename columns
2130 table.rename_column("id", "objectId")
2131 table.rename_column("parent", "parentObjectId")
2132 table.rename_column("merge_peak_sky", "sky_object")
2134 # Add tract and patch columns
2135 table["tract"] = ref.dataId["tract"]
2136 table["patch"] = ref.dataId["patch"]
2138 tables.append(table)
2139 outputTable = astropy.table.vstack(tables, join_type="exact")
2140 butlerQC.put(pipeBase.Struct(outputCatalog=outputTable), outputRefs)