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