Coverage for python/lsst/pipe/tasks/postprocess.py: 57%

901 statements  

« prev     ^ index     » next       coverage.py v7.16.2, created at 2026-10-16 09:04 +0000

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: 533 ↛ 534line 533 didn't jump to line 534 because the condition on line 533 was never true

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): 550 ↛ 553line 550 didn't jump to line 553 because the condition on line 550 was always true

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: 566 ↛ 568line 566 didn't jump to line 568 because the condition on line 566 was always true

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): 572 ↛ 573line 572 didn't jump to line 573 because the condition on line 572 was never true

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: 710 ↛ 716line 710 didn't jump to line 716 because the condition on line 710 was always true

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: 770 ↛ 771line 770 didn't jump to line 771 because the condition on line 770 was never true

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: 778 ↛ 789line 778 didn't jump to line 789 because the condition on line 778 was always true

779 for key in self.config.columnsFromDataId: 

780 if key in dataId: 780 ↛ 787line 780 didn't jump to line 787 because the condition on line 780 was always true

781 if key == "detector": 781 ↛ 783line 781 didn't jump to line 783 because the condition on line 781 was never true

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: 789 ↛ 794line 789 didn't jump to line 794 because the condition on line 789 was always true

790 if df.index.name != self.config.primaryKey and self.config.primaryKey in df: 790 ↛ 791line 790 didn't jump to line 791 because the condition on line 790 was never true

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): 948 ↛ 949line 948 didn't jump to line 949 because the condition on line 948 was never true

949 funcDict_in = funcs.funcDict 

950 elif isinstance(funcs, dict): 950 ↛ 952line 950 didn't jump to line 952 because the condition on line 950 was always true

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: 960 ↛ 961line 960 didn't jump to line 961 because the condition on line 960 was never true

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): 982 ↛ 983line 982 didn't jump to line 983 because the condition on line 982 was never true

983 if band not in outputBands: 

984 continue 

985 # This is something like a ReferenceBand that has configurable bands 

986 elif hasattr(func, "bands"): 986 ↛ 991line 986 didn't jump to line 991 because the condition on line 986 was never true

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"): 992 ↛ 993line 992 didn't jump to line 993 because the condition on line 992 was never true

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): 1029 ↛ 1030line 1029 didn't jump to line 1030 because the condition on line 1029 was never true

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: 1045 ↛ 1046line 1045 didn't jump to line 1046 because the condition on line 1045 was never true

1046 noDupCols.remove(self.config.primaryKey) 

1047 if dataId and self.config.columnsFromDataId: 1047 ↛ 1049line 1047 didn't jump to line 1049 because the condition on line 1047 was always true

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): 1059 ↛ 1072line 1059 didn't jump to line 1072 because the condition on line 1059 was always true

1060 # Allow astropy table inputs to already have the output index 

1061 if name_index not in df_dataset.colnames: 1061 ↛ 1062line 1061 didn't jump to line 1062 because the condition on line 1061 was never true

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: 1078 ↛ 1079line 1078 didn't jump to line 1079 because the condition on line 1078 was never true

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, [])}: 1081 ↛ 1082line 1081 didn't jump to line 1082 because the loop on line 1081 never started

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: 1094 ↛ 1101line 1094 didn't jump to line 1101 because the condition on line 1094 was always true

1095 columns_drop = [column for column in self.config.columnsFromDataId if column in result.df] 

1096 if columns_drop: 1096 ↛ 1101line 1096 didn't jump to line 1101 because the condition on line 1096 was always true

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 default=0.25, 

1247 ) 

1248 min_n_good_to_shift_flux_ratio = pexConfig.Field[int]( 

1249 doc="Minimum number of objects with PSF S/N > good_sn_min and with " 

1250 " size larger than max_reff_compact to use to compute the median " 

1251 " PSF-to-model flux ratio, which is assumed to be 1 otherwise." 

1252 " If this value is not >0, the median flux ratio will be kept 1.", 

1253 default=0, 

1254 ) 

1255 

1256 def __call__(self, table: astropy.table.Table) -> dict[str, NDArray]: 

1257 size_column = self.size_column 

1258 size_model = np.sqrt( 

1259 0.5*(table[size_column.format(axis='x')]**2 + table[size_column.format(axis='y')]**2) 

1260 ) 

1261 small = size_model < self.max_reff_compact 

1262 n_obj = len(table) 

1263 band_mappings_to_process = {band: [band] for band in self.bands} 

1264 band_mappings_to_process.update({k: v.split(",") for k, v in self.bands_combined.items()}) 

1265 output = {} 

1266 for output_band, input_bands in band_mappings_to_process.items(): 

1267 if len(input_bands) > 1: 

1268 flux_psf, fluxerr_psf_sq, flux_model, fluxerr_model_sq = ( 

1269 np.zeros(n_obj, dtype=float) for _ in range(4)) 

1270 for input_band in input_bands: 

1271 flux_psf_b, fluxerr_psf_b, flux_model_b, fluxerr_model_b = self._get_fluxes( 

1272 table, band=input_band) 

1273 # There's no point adding S/N < 0 fluxes 

1274 good = np.isfinite(flux_psf_b) & np.isfinite(flux_model_b) & ( 

1275 flux_psf_b > 0) & (flux_model_b > 0) & (fluxerr_psf_b > 0) & (fluxerr_model_b > 0) 

1276 flux_psf[good] += flux_psf_b[good] 

1277 flux_model[good] += flux_model_b[good] 

1278 fluxerr_psf_sq[good] += fluxerr_psf_b[good]**2 

1279 fluxerr_model_sq[good] += fluxerr_model_b[good]**2 

1280 fluxerr_psf = np.sqrt(fluxerr_psf_sq) 

1281 fluxerr_model = np.sqrt(fluxerr_model_sq) 

1282 fluxerr_psf[fluxerr_psf == 0] = np.inf 

1283 fluxerr_model[fluxerr_model == 0] = np.inf 

1284 else: 

1285 flux_psf, fluxerr_psf, flux_model, fluxerr_model = self._get_fluxes( 

1286 table, band=input_bands[0]) 

1287 

1288 psf_sn = flux_psf/fluxerr_psf 

1289 flux_ratio = np.array(flux_psf / flux_model) 

1290 

1291 if self.min_n_good_to_shift_flux_ratio > 0: 1291 ↛ 1298line 1291 didn't jump to line 1298 because the condition on line 1291 was always true

1292 good = small & (psf_sn > self.good_sn_min) 

1293 # Attempt to correct any flux-independent systematic offset 

1294 # Might need to be a function of S/N 

1295 if np.sum(good == True) > self.min_n_good_to_shift_flux_ratio: # noqa: E712 1295 ↛ 1298line 1295 didn't jump to line 1298 because the condition on line 1295 was always true

1296 flux_ratio *= 1./np.nanmedian(flux_ratio[good]) 

1297 

1298 flux_ratio_err = np.sqrt( 

1299 (fluxerr_psf/flux_model)**2 + (fluxerr_model*fluxerr_psf/flux_model**2)**2 

1300 ) 

1301 extendedness = (1 - flux_ratio) + self.fluxerr_coefficent*flux_ratio_err 

1302 extendedness *= np.sqrt(size_model/self.max_reff_compact) 

1303 extendedness[(extendedness < 0) & (extendedness > -np.inf)] = 0 

1304 # Make it sigmoid-like with a stretch 

1305 stretch = self.fluxerr_stretch 

1306 extendedness *= stretch 

1307 extendedness = np.clip((stretch + 1)/stretch*extendedness/(1 + extendedness), 0, 1) 

1308 

1309 column_out = self.output_column.format(band=output_band) 

1310 output[column_out] = extendedness 

1311 return output 

1312 

1313 

1314class ConsolidateObjectTableConnections(pipeBase.PipelineTaskConnections, 

1315 dimensions=("tract", "skymap")): 

1316 inputCatalogs = connectionTypes.Input( 

1317 doc="Per-Patch objectTables conforming to the standard data model.", 

1318 name="objectTable", 

1319 storageClass="ArrowAstropy", 

1320 dimensions=("tract", "patch", "skymap"), 

1321 multiple=True, 

1322 deferLoad=True, 

1323 ) 

1324 outputCatalog = connectionTypes.Output( 

1325 doc="Pre-tract horizontal concatenation of the input objectTables", 

1326 name="objectTable_tract", 

1327 storageClass="ArrowAstropy", 

1328 dimensions=("tract", "skymap"), 

1329 ) 

1330 

1331 

1332class ConsolidateObjectTableConfig(pipeBase.PipelineTaskConfig, 

1333 pipelineConnections=ConsolidateObjectTableConnections): 

1334 actions = ConfigurableActionStructField[ComputeColumnsAction]( 

1335 doc="Actions to add columns to the final object table", 

1336 ) 

1337 coaddName = pexConfig.Field[str]( 

1338 default="deep", 

1339 doc="Name of coadd" 

1340 ) 

1341 

1342 def setDefaults(self): 

1343 super().setDefaults() 

1344 self.actions.extendedness = ModelExtendednessColumnAction() 

1345 

1346 

1347class ConsolidateObjectTableTask(pipeBase.PipelineTask): 

1348 """Write patch-merged source tables to a tract-level DataFrame Parquet file. 

1349 

1350 Concatenates `objectTable` list into a per-visit `objectTable_tract`. 

1351 """ 

1352 _DefaultName = "consolidateObjectTable" 

1353 ConfigClass = ConsolidateObjectTableConfig 

1354 

1355 inputDataset = "objectTable" 

1356 outputDataset = "objectTable_tract" 

1357 

1358 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

1359 inputs = butlerQC.get(inputRefs) 

1360 self.log.info("Concatenating %s per-patch Object Tables", 

1361 len(inputs["inputCatalogs"])) 

1362 table = TableVStack.vstack_handles(inputs["inputCatalogs"]) 

1363 for action in self.config.actions: 

1364 computed = action(table) 

1365 for key, values in computed.items(): 

1366 table[key] = values.astype(np.float32) if values.dtype == np.float64 else values 

1367 butlerQC.put(pipeBase.Struct(outputCatalog=table), outputRefs) 

1368 

1369 

1370class TransformSourceTableConnections(pipeBase.PipelineTaskConnections, 

1371 defaultTemplates={"catalogType": ""}, 

1372 dimensions=("instrument", "visit", "detector")): 

1373 

1374 inputCatalog = connectionTypes.Input( 

1375 doc="Wide input catalog of sources produced by WriteSourceTableTask", 

1376 name="{catalogType}source", 

1377 storageClass="DataFrame", 

1378 dimensions=("instrument", "visit", "detector"), 

1379 deferLoad=True 

1380 ) 

1381 outputCatalog = connectionTypes.Output( 

1382 doc="Narrower, per-detector Source Table transformed and converted per a " 

1383 "specified set of functors", 

1384 name="{catalogType}sourceTable", 

1385 storageClass="ArrowAstropy", 

1386 dimensions=("instrument", "visit", "detector") 

1387 ) 

1388 

1389 

1390class TransformSourceTableConfig(TransformCatalogBaseConfig, 

1391 pipelineConnections=TransformSourceTableConnections): 

1392 

1393 def setDefaults(self): 

1394 super().setDefaults() 

1395 self.functorFile = os.path.join("$PIPE_TASKS_DIR", "schemas", "Source.yaml") 

1396 self.primaryKey = "sourceId" 

1397 self.columnsFromDataId = ["visit", "detector", "band", "physical_filter"] 

1398 

1399 

1400class TransformSourceTableTask(TransformCatalogBaseTask): 

1401 """Transform/standardize a source catalog 

1402 """ 

1403 _DefaultName = "transformSourceTable" 

1404 ConfigClass = TransformSourceTableConfig 

1405 

1406 

1407class ConsolidateVisitSummaryConnections(pipeBase.PipelineTaskConnections, 

1408 dimensions=("instrument", "visit",), 

1409 defaultTemplates={"calexpType": ""}): 

1410 camera = connectionTypes.PrerequisiteInput( 

1411 doc="Camera geometry.", 

1412 name="camera", 

1413 dimensions=("instrument",), 

1414 storageClass="Camera", 

1415 isCalibration=True, 

1416 ) 

1417 calexp = connectionTypes.Input( 

1418 doc="Processed exposures used for metadata. These are legacy Exposures " 

1419 "or lsst.images.VisitImages depending on input_image_type.", 

1420 name="calexp", 

1421 storageClass="ExposureF", 

1422 dimensions=("instrument", "visit", "detector"), 

1423 deferLoad=True, 

1424 multiple=True, 

1425 ) 

1426 visitSummary = connectionTypes.Output( 

1427 doc=("Per-visit consolidated exposure metadata. These catalogs use " 

1428 "detector id for the id and are sorted for fast lookups of a " 

1429 "detector."), 

1430 name="visitSummary", 

1431 storageClass="ExposureCatalog", 

1432 dimensions=("instrument", "visit"), 

1433 ) 

1434 visitSummarySchema = connectionTypes.InitOutput( 

1435 doc="Schema of the visitSummary catalog", 

1436 name="visitSummary_schema", 

1437 storageClass="ExposureCatalog", 

1438 ) 

1439 # RefitPointingTask options and connections use snake_case for consistency 

1440 # with options of the same name in UpdateVisitSummaryTask, which seems more 

1441 # important than consistently using camelCase within this task. 

1442 visit_geometry = connectionTypes.Output( 

1443 doc="Updated visit[, detector] regions that can be used to update butler dimensions records.", 

1444 name="visit_geometry", 

1445 dimensions=("instrument", "visit"), 

1446 storageClass="VisitGeometry", 

1447 ) 

1448 

1449 def __init__(self, *, config=None): 

1450 if self.config.do_refit_pointing: 

1451 self.camera = dataclasses.replace(self.camera, dimensions=config.cameraDimensions) 

1452 else: 

1453 del self.camera 

1454 if not self.config.do_write_visit_geometry: 

1455 del self.visit_geometry 

1456 if self.config.input_image_type == "future": 

1457 self.calexp = dataclasses.replace(self.calexp, storageClass="VisitImage") 

1458 

1459 

1460class ConsolidateVisitSummaryConfig(pipeBase.PipelineTaskConfig, 

1461 pipelineConnections=ConsolidateVisitSummaryConnections): 

1462 """Config for ConsolidateVisitSummaryTask""" 

1463 

1464 full = pexConfig.Field( 

1465 "Whether to propate all exposure components. " 

1466 "This adds PSF, aperture correction map, transmission curve, and detector, which can increase file " 

1467 "size by more than factor of 10, but it makes the visit summaries produced by this task fully usable" 

1468 "by tasks that were designed to run downstream of lsst.drp.tasks.UpdateVisitSummaryTask.", 

1469 dtype=bool, 

1470 default=False, 

1471 ) 

1472 refitPointing = pexConfig.ConfigurableField( 

1473 "A subtask for refitting the boresight pointing and orientation, " 

1474 "and using those to produce new regions for butler dimensions.", 

1475 target=RefitPointingTask, 

1476 ) 

1477 cameraDimensions = pexConfig.ListField( 

1478 "The dimensions of the 'camera' prerequisite input connection.", 

1479 dtype=str, 

1480 default=["instrument"], 

1481 ) 

1482 do_refit_pointing = pexConfig.Field( 

1483 "Whether to re-fit the pointing model.", 

1484 dtype=bool, 

1485 default=True, 

1486 ) 

1487 do_write_visit_geometry = pexConfig.Field( 

1488 "Whether to write refit-pointing regions that can be used to update butler dimension records.", 

1489 dtype=bool, 

1490 default=True, 

1491 ) 

1492 input_image_type = pexConfig.ChoiceField[str]( 

1493 "Which image type to read the per-detector metadata from.", 

1494 allowed={ 

1495 "legacy": "Read components of an lsst.afw.image.Exposure.", 

1496 "future": "Read components of an lsst.images.VisitImage.", 

1497 }, 

1498 optional=False, 

1499 default="legacy", 

1500 ) 

1501 

1502 def validate(self): 

1503 super().validate() 

1504 if self.do_write_visit_geometry and not self.do_refit_pointing: 1504 ↛ 1505line 1504 didn't jump to line 1505 because the condition on line 1504 was never true

1505 raise ValueError("Cannot write visit_geometry without refitting the pointing.") 

1506 if self.full and self.input_image_type == "future": 

1507 raise ValueError("full=True requires input_image_type='legacy', because an " 

1508 "lsst.images.VisitImage has no transmission curve.") 

1509 

1510 def setDefaults(self): 

1511 super().setDefaults() 

1512 # This prevents a conflict with the fields added by the 

1513 # UpdateVisitSummaryTask invocation of RefitPointingTask. 

1514 self.refitPointing.schema_prefix = "preliminary_" 

1515 

1516 

1517class ConsolidateVisitSummaryTask(pipeBase.PipelineTask): 

1518 """Task to consolidate per-detector visit metadata. 

1519 

1520 This task aggregates the following metadata from all the detectors in a 

1521 single visit into an exposure catalog: 

1522 - The visitInfo. 

1523 - The wcs. 

1524 - The photoCalib. 

1525 - The physical_filter and band (if available). 

1526 - The PSF model. 

1527 - The aperture correction map. 

1528 - The transmission curve. 

1529 - The psf size, shape, and effective area at the center of the detector. 

1530 - The corners of the bounding box in right ascension/declination. 

1531 

1532 Tests for this task are performed in ci_hsc_gen3. 

1533 """ 

1534 _DefaultName = "consolidateVisitSummary" 

1535 ConfigClass = ConsolidateVisitSummaryConfig 

1536 

1537 def __init__(self, **kwargs): 

1538 super().__init__(**kwargs) 

1539 self.schema = afwTable.ExposureTable.makeMinimalSchema() 

1540 self.schema.addField("visit", type="L", doc="Visit number") 

1541 self.schema.addField("physical_filter", type="String", size=32, doc="Physical filter") 

1542 self.schema.addField("band", type="String", size=32, doc="Name of band") 

1543 ExposureSummaryStats.update_schema(self.schema) 

1544 if self.config.do_refit_pointing: 

1545 self.makeSubtask("refitPointing", schema=self.schema) 

1546 self.visitSummarySchema = afwTable.ExposureCatalog(self.schema) 

1547 

1548 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

1549 handles = butlerQC.get(inputRefs.calexp) 

1550 visit = handles[0].dataId["visit"] 

1551 

1552 self.log.debug("Concatenating metadata from %d per-detector calexps (visit %d)", 

1553 len(handles), visit) 

1554 

1555 camera = butlerQC.get(inputRefs.camera) if self.config.do_refit_pointing else None 

1556 

1557 result = pipeBase.Struct() 

1558 try: 

1559 self.run(visit=visit, handles=handles, camera=camera, result=result) 

1560 except pipeBase.AlgorithmError as e: 

1561 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

1562 e, self, result.visitSummary, log=self.log 

1563 ) 

1564 butlerQC.put(result, outputRefs) 

1565 raise error from e 

1566 

1567 butlerQC.put(result, outputRefs) 

1568 

1569 def run(self, *, visit, handles, camera=None, result=None): 

1570 """Make a combined exposure catalog from a list of handles. 

1571 These handles must point to exposures with wcs, summaryStats, 

1572 and other visit metadata. 

1573 

1574 Parameters 

1575 ---------- 

1576 visit : `int` 

1577 Visit identification number. 

1578 handles : `list` of `lsst.daf.butler.DeferredDatasetHandle` 

1579 List of handles in visit, to `lsst.afw.image.Exposure` or 

1580 `lsst.images.VisitImage` datasets depending on 

1581 ``config.input_image_type``. 

1582 camera : `lsst.afw.cameraGeom.Camera`, optional 

1583 Camera geometry. Required if and only if 

1584 ``do_refit_pointing=True``. 

1585 result : `lsst.pipe.base.Struct`, optional 

1586 Output struct to modify in-place. 

1587 

1588 Returns 

1589 ------- 

1590 result : `lsst.pipe.base.Struct` 

1591 Struct with the following attributes: 

1592 

1593 - ``visitSummary`` (`lsst.afw.table.ExposureCatalog`): an Exposure 

1594 catalog with per-detector summary information. 

1595 - ``visit_geometry`` (`lsst.obs.base.visit_geometry.VisitGeometry`): 

1596 Regions that can be used to update butler dimension regions for 

1597 this visit. Only present if ``do_refit_pointing=True``. 

1598 """ 

1599 cat = afwTable.ExposureCatalog(self.schema) 

1600 cat.resize(len(handles)) 

1601 

1602 cat["visit"] = visit 

1603 

1604 # These are common for all detectors and should not 

1605 # be read in more than once. 

1606 visitInfo = None 

1607 filterLabel = None 

1608 

1609 for i, dataRef in enumerate(handles): 

1610 if self.config.input_image_type == "future": 

1611 components = self._readFutureComponents(dataRef, readVisitInfo=visitInfo is None) 

1612 else: 

1613 components = self._readLegacyComponents(dataRef, readVisitInfo=visitInfo is None) 

1614 if visitInfo is None: 1614 ↛ 1617line 1614 didn't jump to line 1617 because the condition on line 1614 was always true

1615 visitInfo = components.visitInfo 

1616 filterLabel = components.filterLabel 

1617 summaryStats = components.summaryStats 

1618 

1619 rec = cat[i] 

1620 rec.setBBox(components.bbox) 

1621 rec.setVisitInfo(visitInfo) 

1622 rec.setWcs(components.wcs) 

1623 rec.setPhotoCalib(components.photoCalib) 

1624 rec.setValidPolygon(components.validPolygon) 

1625 

1626 if self.config.full: 1626 ↛ 1627line 1626 didn't jump to line 1627 because the condition on line 1626 was never true

1627 rec.setPsf(dataRef.get(component="psf")) 

1628 rec.setApCorrMap(dataRef.get(component="apCorrMap")) 

1629 rec.setTransmissionCurve(dataRef.get(component="transmissionCurve")) 

1630 

1631 rec["physical_filter"] = filterLabel.physicalLabel if filterLabel.hasPhysicalLabel() else "" 

1632 rec["band"] = filterLabel.bandLabel if filterLabel.hasBandLabel() else "" 

1633 rec.setId(dataRef.dataId["detector"]) 

1634 summaryStats.update_record(rec) 

1635 

1636 if not cat: 1636 ↛ 1637line 1636 didn't jump to line 1637 because the condition on line 1636 was never true

1637 raise pipeBase.NoWorkFound( 

1638 "No detectors had sufficient information to make a visit summary row." 

1639 ) 

1640 

1641 metadata = dafBase.PropertyList() 

1642 metadata.add("COMMENT", "Catalog id is detector id, sorted.") 

1643 # We are looping over existing handles, so the following is true 

1644 metadata.add("COMMENT", "Only detectors with data have entries.") 

1645 cat.setMetadata(metadata) 

1646 

1647 cat.sort() 

1648 

1649 if result is None: 1649 ↛ 1651line 1649 didn't jump to line 1651 because the condition on line 1649 was always true

1650 result = pipeBase.Struct() 

1651 result.visitSummary = cat 

1652 

1653 if self.config.do_refit_pointing: 1653 ↛ 1654line 1653 didn't jump to line 1654 because the condition on line 1653 was never true

1654 refitPointingResult = self.refitPointing.run(catalog=cat, camera=camera) 

1655 result.visit_geometry = refitPointingResult.regions 

1656 

1657 return result 

1658 

1659 @staticmethod 

1660 def _readLegacyComponents(handle, readVisitInfo): 

1661 """Read the metadata of one detector from an 

1662 `lsst.afw.image.Exposure`. 

1663 

1664 Parameters 

1665 ---------- 

1666 handle : `lsst.daf.butler.DeferredDatasetHandle` 

1667 Handle to the exposure. 

1668 readVisitInfo : `bool` 

1669 Whether to read the visit info and filter label, which are the same 

1670 for every detector in a visit. 

1671 

1672 Returns 

1673 ------- 

1674 components : `lsst.pipe.base.Struct` 

1675 Struct with ``visitInfo``, ``filterLabel`` (`None` unless 

1676 ``readVisitInfo``), ``summaryStats``, ``wcs``, ``photoCalib``, 

1677 ``bbox`` and ``validPolygon`` attributes. 

1678 """ 

1679 return pipeBase.Struct( 

1680 visitInfo=handle.get(component="visitInfo") if readVisitInfo else None, 

1681 filterLabel=handle.get(component="filter") if readVisitInfo else None, 

1682 summaryStats=handle.get(component="summaryStats"), 

1683 wcs=handle.get(component="wcs"), 

1684 photoCalib=handle.get(component="photoCalib"), 

1685 bbox=handle.get(component="bbox"), 

1686 validPolygon=handle.get(component="validPolygon"), 

1687 ) 

1688 

1689 @staticmethod 

1690 def _readFutureComponents(handle, readVisitInfo): 

1691 """Read the metadata of one detector from an `lsst.images.VisitImage` 

1692 and convert it to legacy types. 

1693 

1694 Parameters 

1695 ---------- 

1696 handle : `lsst.daf.butler.DeferredDatasetHandle` 

1697 Handle to the visit image. 

1698 readVisitInfo : `bool` 

1699 Whether to read the visit info and filter label, which are the same 

1700 for every detector in a visit. 

1701 

1702 Returns 

1703 ------- 

1704 components : `lsst.pipe.base.Struct` 

1705 Struct with the same attributes as `_readLegacyComponents`. 

1706 """ 

1707 names = ["summary_stats", "sky_projection", "photometric_scaling", "unit", "bbox", "bounds"] 

1708 if readVisitInfo: 1708 ↛ 1711line 1708 didn't jump to line 1711 because the condition on line 1708 was always true

1709 names.append("obs_info") 

1710 # One multi-component read opens the file once and skips the pixels. 

1711 read = handle.get(component="components", parameters={"components": names}) 

1712 if read["photometric_scaling"] is not None: 

1713 photoCalib = read["photometric_scaling"].to_legacy_photo_calib(read["unit"]) 

1714 else: 

1715 photoCalib = BaseField.make_legacy_photo_calib(read["unit"]) 

1716 visitInfo = None 

1717 filterLabel = None 

1718 if readVisitInfo: 1718 ↛ 1721line 1718 didn't jump to line 1721 because the condition on line 1718 was always true

1719 visitInfo = MakeRawVisitInfoViaObsInfo.observationInfo2visitInfo(read["obs_info"]) 

1720 filterLabel = FilterLabel(band=handle.dataId["band"], physical=handle.dataId["physical_filter"]) 

1721 bounds = read["bounds"] 

1722 return pipeBase.Struct( 

1723 visitInfo=visitInfo, 

1724 filterLabel=filterLabel, 

1725 summaryStats=read["summary_stats"].to_legacy(), 

1726 wcs=read["sky_projection"].to_legacy(), 

1727 photoCalib=photoCalib, 

1728 bbox=read["bbox"].to_legacy(), 

1729 validPolygon=bounds.to_legacy() if isinstance(bounds, Polygon) else None, 

1730 ) 

1731 

1732 

1733class ConsolidateSourceTableConnections(pipeBase.PipelineTaskConnections, 

1734 defaultTemplates={"catalogType": ""}, 

1735 dimensions=("instrument", "visit")): 

1736 inputCatalogs = connectionTypes.Input( 

1737 doc="Input per-detector Source Tables", 

1738 name="{catalogType}sourceTable", 

1739 storageClass="ArrowAstropy", 

1740 dimensions=("instrument", "visit", "detector"), 

1741 multiple=True, 

1742 deferLoad=True, 

1743 ) 

1744 outputCatalog = connectionTypes.Output( 

1745 doc="Per-visit concatenation of Source Table", 

1746 name="{catalogType}sourceTable_visit", 

1747 storageClass="ArrowAstropy", 

1748 dimensions=("instrument", "visit") 

1749 ) 

1750 

1751 

1752class ConsolidateSourceTableConfig(pipeBase.PipelineTaskConfig, 

1753 pipelineConnections=ConsolidateSourceTableConnections): 

1754 pass 

1755 

1756 

1757class ConsolidateSourceTableTask(pipeBase.PipelineTask): 

1758 """Concatenate `sourceTable` list into a per-visit `sourceTable_visit` 

1759 """ 

1760 _DefaultName = "consolidateSourceTable" 

1761 ConfigClass = ConsolidateSourceTableConfig 

1762 

1763 inputDataset = "sourceTable" 

1764 outputDataset = "sourceTable_visit" 

1765 

1766 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

1767 # Docstring inherited. 

1768 detectorOrder = [ref.dataId["detector"] for ref in inputRefs.inputCatalogs] 

1769 detectorOrder.sort() 

1770 inputRefs = reorderRefs(inputRefs, detectorOrder, dataIdKey="detector") 

1771 inputs = butlerQC.get(inputRefs) 

1772 self.log.info("Concatenating %s per-detector Source Tables", 

1773 len(inputs["inputCatalogs"])) 

1774 table = TableVStack.vstack_handles(inputs["inputCatalogs"]) 

1775 butlerQC.put(pipeBase.Struct(outputCatalog=table), outputRefs) 

1776 

1777 

1778class MakeCcdVisitTableConnections(pipeBase.PipelineTaskConnections, 

1779 dimensions=("instrument",), 

1780 defaultTemplates={"calexpType": ""}): 

1781 visitSummaryRefs = connectionTypes.Input( 

1782 doc="Data references for per-visit consolidated exposure metadata", 

1783 name="finalVisitSummary", 

1784 storageClass="ExposureCatalog", 

1785 dimensions=("instrument", "visit"), 

1786 multiple=True, 

1787 deferLoad=True, 

1788 ) 

1789 outputCatalog = connectionTypes.Output( 

1790 doc="CCD and Visit metadata table", 

1791 name="ccdVisitTable", 

1792 storageClass="ArrowAstropy", 

1793 dimensions=("instrument",) 

1794 ) 

1795 

1796 

1797class MakeCcdVisitTableConfig(pipeBase.PipelineTaskConfig, 

1798 pipelineConnections=MakeCcdVisitTableConnections): 

1799 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

1800 

1801 

1802class MakeCcdVisitTableTask(pipeBase.PipelineTask): 

1803 """Produce a `ccdVisitTable` from the visit summary exposure catalogs. 

1804 """ 

1805 _DefaultName = "makeCcdVisitTable" 

1806 ConfigClass = MakeCcdVisitTableConfig 

1807 

1808 def run(self, visitSummaryRefs): 

1809 """Make a table of ccd information from the visit summary catalogs. 

1810 

1811 Parameters 

1812 ---------- 

1813 visitSummaryRefs : `list` of `lsst.daf.butler.DeferredDatasetHandle` 

1814 List of DeferredDatasetHandles pointing to exposure catalogs with 

1815 per-detector summary information. 

1816 

1817 Returns 

1818 ------- 

1819 result : `~lsst.pipe.base.Struct` 

1820 Results struct with attribute: 

1821 

1822 ``outputCatalog`` 

1823 Catalog of ccd and visit information. 

1824 """ 

1825 ccdEntries = [] 

1826 for visitSummaryRef in visitSummaryRefs: 

1827 visitSummary = visitSummaryRef.get() 

1828 if not visitSummary: 

1829 continue 

1830 visitInfo = visitSummary[0].getVisitInfo() 

1831 

1832 # Strip provenance to prevent merge confusion. 

1833 strip_provenance_from_fits_header(visitSummary.metadata) 

1834 

1835 ccdEntry = {} 

1836 summaryTable = visitSummary.asAstropy() 

1837 selectColumns = ["id", "visit", "physical_filter", "band", "ra", "dec", 

1838 "pixelScale", "zenithDistance", 

1839 "expTime", "zeroPoint", "psfSigma", 

1840 "psfArea", "psfIxx", "psfIyy", "psfIxy", 

1841 "skyBg", "skyNoise", 

1842 "psfAdaptiveThresholdValue", "psfAdaptiveIncludeThresholdMultiplier", 

1843 "nShapeletsStar", "shapeletsOnlyIqScore", "shapeletsIqScore", 

1844 "centroidDiffShapeletsVsSlotMedian", 

1845 "shapeletsStarEMedian", "shapeletsStarUnNormalizedEMedian", 

1846 "refCatSourceDensity", "astromOffsetMean", "astromOffsetStd", 

1847 "nPsfStar", "psfStarDeltaE1Median", "psfStarDeltaE2Median", 

1848 "psfStarDeltaE1Scatter", "psfStarDeltaE2Scatter", 

1849 "psfStarDeltaSizeMedian", "psfStarDeltaSizeScatter", 

1850 "psfStarScaledDeltaSizeScatter", "psfTraceRadiusDelta", 

1851 "psfApFluxDelta", "psfApCorrSigmaScaledDelta", 

1852 "maxDistToNearestPsf", "starEMedian", "starUnNormalizedEMedian", 

1853 "starComa1Median", "starComa2Median", 

1854 "starTrefoil1Median", "starTrefoil2Median", 

1855 "starKurtosisMedian", "starE41Median", "starE42Median", 

1856 "effTime", "effTimePsfSigmaScale", 

1857 "effTimeSkyBgScale", "effTimeZeroPointScale", 

1858 "magLim"] 

1859 ccdEntry = summaryTable[selectColumns] 

1860 # 'visit' is the human readable visit number. 

1861 # 'visitId' is the key to the visitId table. They are the same. 

1862 # Technically you should join to get the visit from the visit 

1863 # table. 

1864 ccdEntry.rename_column("visit", "visitId") 

1865 ccdEntry.rename_column("id", "detectorId") 

1866 

1867 # RFC-924: Temporarily keep a duplicate "decl" entry for backwards 

1868 # compatibility. To be removed after September 2023. 

1869 ccdEntry["decl"] = ccdEntry["dec"] 

1870 

1871 ccdEntry["ccdVisitId"] = [ 

1872 self.config.idGenerator.apply( 

1873 visitSummaryRef.dataId, 

1874 detector=detector_id, 

1875 is_exposure=False, 

1876 ).catalog_id # The "catalog ID" here is the ccdVisit ID 

1877 # because it's usually the ID for a whole catalog 

1878 # with a {visit, detector}, and that's the main 

1879 # use case for IdGenerator. This usage for a 

1880 # summary table is rare. 

1881 for detector_id in summaryTable["id"] 

1882 ] 

1883 ccdEntry["detector"] = summaryTable["id"] 

1884 ccdEntry["seeing"] = ( 

1885 visitSummary["psfSigma"] * visitSummary["pixelScale"] * np.sqrt(8 * np.log(2)) 

1886 ) 

1887 ccdEntry["skyRotation"] = visitInfo.getBoresightRotAngle().asDegrees() 

1888 ccdEntry["expMidpt"] = np.datetime64(visitInfo.date.nsecs(scale=dafBase.DateTime.TAI), "ns") 

1889 ccdEntry["expMidptMJD"] = visitInfo.getDate().get(dafBase.DateTime.MJD) 

1890 expTime = visitInfo.getExposureTime() 

1891 ccdEntry["obsStart"] = ( 

1892 ccdEntry["expMidpt"] - 0.5 * np.timedelta64(int(expTime * 1E9), "ns") 

1893 ) 

1894 expTime_days = expTime / (60*60*24) 

1895 ccdEntry["obsStartMJD"] = ccdEntry["expMidptMJD"] - 0.5 * expTime_days 

1896 ccdEntry["darkTime"] = visitInfo.getDarkTime() 

1897 ccdEntry["xSize"] = summaryTable["bbox_max_x"] - summaryTable["bbox_min_x"] 

1898 ccdEntry["ySize"] = summaryTable["bbox_max_y"] - summaryTable["bbox_min_y"] 

1899 ccdEntry["llcra"] = summaryTable["raCorners"][:, 0] 

1900 ccdEntry["llcdec"] = summaryTable["decCorners"][:, 0] 

1901 ccdEntry["ulcra"] = summaryTable["raCorners"][:, 1] 

1902 ccdEntry["ulcdec"] = summaryTable["decCorners"][:, 1] 

1903 ccdEntry["urcra"] = summaryTable["raCorners"][:, 2] 

1904 ccdEntry["urcdec"] = summaryTable["decCorners"][:, 2] 

1905 ccdEntry["lrcra"] = summaryTable["raCorners"][:, 3] 

1906 ccdEntry["lrcdec"] = summaryTable["decCorners"][:, 3] 

1907 # These columns are only added on certain updateVisitSummary 

1908 # configurations, and we want to rename them to camelCase and 

1909 # convert to arcseconds here. 

1910 for inName, outName in { 

1911 "wcs_corner_max_offset": "wcsCornerMaxOffset", 

1912 "wcs_detector_pointing_residual": "wcsDetectorPointingResidual", 

1913 "wcs_visit_pointing_residual": "wcsVisitPointingResidual", 

1914 "preliminary_wcs_detector_pointing_residual": "wcsPreliminaryDetectorPointingResidual", 

1915 "preliminary_wcs_visit_pointing_residual": "wcsPreliminaryVisitPointingResidual", 

1916 "wcs_sip_delta_sky": "wcsSIPDeltaSky", 

1917 }.items(): 

1918 if inName in summaryTable.columns: 

1919 inCol = summaryTable[inName] 

1920 ccdEntry[outName] = astropy.table.Column( 

1921 (np.asarray(inCol) * inCol.unit).to_value(astropy.units.arcsec), 

1922 unit=astropy.units.arcsec 

1923 ) 

1924 if "wcs_sip_delta_pixel" in summaryTable.columns: 

1925 ccdEntry["wcsSIPDeltaPixel"] = summaryTable["wcs_sip_delta_pixel"] 

1926 

1927 # TODO: DM-30618, Add raftName, nExposures, ccdTemp, binX, binY, 

1928 # and flags, and decide if WCS, and llcx, llcy, ulcx, ulcy, etc. 

1929 # values are actually wanted. 

1930 ccdEntries.append(ccdEntry) 

1931 

1932 outputCatalog = astropy.table.vstack(ccdEntries, join_type="exact") 

1933 return pipeBase.Struct(outputCatalog=outputCatalog) 

1934 

1935 

1936class MakeVisitTableConnections(pipeBase.PipelineTaskConnections, 

1937 dimensions=("instrument",), 

1938 defaultTemplates={"calexpType": ""}): 

1939 visitSummaries = connectionTypes.Input( 

1940 doc="Per-visit consolidated exposure metadata", 

1941 name="finalVisitSummary", 

1942 storageClass="ExposureCatalog", 

1943 dimensions=("instrument", "visit",), 

1944 multiple=True, 

1945 deferLoad=True, 

1946 ) 

1947 outputCatalog = connectionTypes.Output( 

1948 doc="Visit metadata table", 

1949 name="visitTable", 

1950 storageClass="ArrowAstropy", 

1951 dimensions=("instrument",) 

1952 ) 

1953 

1954 

1955class MakeVisitTableConfig(pipeBase.PipelineTaskConfig, 

1956 pipelineConnections=MakeVisitTableConnections): 

1957 pass 

1958 

1959 

1960class MakeVisitTableTask(pipeBase.PipelineTask): 

1961 """Produce a `visitTable` from the visit summary exposure catalogs. 

1962 """ 

1963 _DefaultName = "makeVisitTable" 

1964 ConfigClass = MakeVisitTableConfig 

1965 

1966 def run(self, visitSummaries): 

1967 """Make a table of visit information from the visit summary catalogs. 

1968 

1969 Parameters 

1970 ---------- 

1971 visitSummaries : `list` of `lsst.afw.table.ExposureCatalog` 

1972 List of exposure catalogs with per-detector summary information. 

1973 Returns 

1974 ------- 

1975 result : `~lsst.pipe.base.Struct` 

1976 Results struct with attribute: 

1977 

1978 ``outputCatalog`` 

1979 Catalog of visit information. 

1980 """ 

1981 visitEntries = [] 

1982 for visitSummary in visitSummaries: 

1983 visitSummary = visitSummary.get() 

1984 if not visitSummary: 

1985 continue 

1986 visitRow = visitSummary[0] 

1987 visitInfo = visitRow.getVisitInfo() 

1988 

1989 visitEntry = {} 

1990 visitEntry["visitId"] = visitRow["visit"] 

1991 visitEntry["visit"] = visitRow["visit"] 

1992 visitEntry["physical_filter"] = visitRow["physical_filter"] 

1993 visitEntry["band"] = visitRow["band"] 

1994 raDec = visitInfo.getBoresightRaDec() 

1995 visitEntry["ra"] = raDec.getRa().asDegrees() 

1996 visitEntry["dec"] = raDec.getDec().asDegrees() 

1997 

1998 # RFC-924: Temporarily keep a duplicate "decl" entry for backwards 

1999 # compatibility. To be removed after September 2023. 

2000 visitEntry["decl"] = visitEntry["dec"] 

2001 

2002 visitEntry["skyRotation"] = visitInfo.getBoresightRotAngle().asDegrees() 

2003 azAlt = visitInfo.getBoresightAzAlt() 

2004 visitEntry["azimuth"] = azAlt.getLongitude().asDegrees() 

2005 visitEntry["altitude"] = azAlt.getLatitude().asDegrees() 

2006 visitEntry["zenithDistance"] = 90 - azAlt.getLatitude().asDegrees() 

2007 visitEntry["airmass"] = visitInfo.getBoresightAirmass() 

2008 expTime = visitInfo.getExposureTime() 

2009 visitEntry["expTime"] = expTime 

2010 visitEntry["expMidpt"] = np.datetime64(visitInfo.date.nsecs(scale=dafBase.DateTime.TAI), "ns") 

2011 visitEntry["expMidptMJD"] = visitInfo.getDate().get(dafBase.DateTime.MJD) 

2012 visitEntry["obsStart"] = visitEntry["expMidpt"] - 0.5 * np.timedelta64(int(expTime * 1E9), "ns") 

2013 expTime_days = expTime / (60*60*24) 

2014 visitEntry["obsStartMJD"] = visitEntry["expMidptMJD"] - 0.5 * expTime_days 

2015 

2016 # Note that these are computed per-visit and stored identically 

2017 # in each row, so we can use the first row values. 

2018 visitEntry["psfTE1e1"] = visitRow["psfTE1e1"] 

2019 visitEntry["psfTE1e2"] = visitRow["psfTE1e2"] 

2020 visitEntry["psfTE1ex"] = visitRow["psfTE1ex"] 

2021 visitEntry["psfTE2e1"] = visitRow["psfTE2e1"] 

2022 visitEntry["psfTE2e2"] = visitRow["psfTE2e2"] 

2023 visitEntry["psfTE2ex"] = visitRow["psfTE2ex"] 

2024 visitEntry["psfTE3e1"] = visitRow["psfTE3e1"] 

2025 visitEntry["psfTE3e2"] = visitRow["psfTE3e2"] 

2026 visitEntry["psfTE3ex"] = visitRow["psfTE3ex"] 

2027 visitEntry["psfTE4e1"] = visitRow["psfTE4e1"] 

2028 visitEntry["psfTE4e2"] = visitRow["psfTE4e2"] 

2029 visitEntry["psfTE4ex"] = visitRow["psfTE4ex"] 

2030 

2031 visitEntries.append(visitEntry) 

2032 

2033 # TODO: DM-30623, Add programId, exposureType, cameraTemp, 

2034 # mirror1Temp, mirror2Temp, mirror3Temp, domeTemp, externalTemp, 

2035 # dimmSeeing, pwvGPS, pwvMW, flags, nExposures. 

2036 

2037 outputCatalog = astropy.table.Table(rows=visitEntries) 

2038 return pipeBase.Struct(outputCatalog=outputCatalog) 

2039 

2040 

2041@deprecated(reason="This task is replaced by lsst.pipe.tasks.ForcedPhotCcdTask. " 

2042 "This task will be removed after v30.", 

2043 version="v29.0", category=FutureWarning) 

2044class WriteForcedSourceTableConnections(pipeBase.PipelineTaskConnections, 

2045 dimensions=("instrument", "visit", "detector", "skymap", "tract")): 

2046 

2047 inputCatalog = connectionTypes.Input( 

2048 doc="Primary per-detector, single-epoch forced-photometry catalog. " 

2049 "By default, it is the output of ForcedPhotCcdTask on calexps", 

2050 name="forced_src", 

2051 storageClass="SourceCatalog", 

2052 dimensions=("instrument", "visit", "detector", "skymap", "tract") 

2053 ) 

2054 inputCatalogDiff = connectionTypes.Input( 

2055 doc="Secondary multi-epoch, per-detector, forced photometry catalog. " 

2056 "By default, it is the output of ForcedPhotCcdTask run on image differences.", 

2057 name="forced_diff", 

2058 storageClass="SourceCatalog", 

2059 dimensions=("instrument", "visit", "detector", "skymap", "tract") 

2060 ) 

2061 outputCatalog = connectionTypes.Output( 

2062 doc="InputCatalogs horizonatally joined on `objectId` in DataFrame parquet format", 

2063 name="mergedForcedSource", 

2064 storageClass="DataFrame", 

2065 dimensions=("instrument", "visit", "detector", "skymap", "tract") 

2066 ) 

2067 

2068 

2069@deprecated(reason="This task is replaced by lsst.pipe.tasks.ForcedPhotCcdTask. " 

2070 "This task will be removed after v30.", 

2071 version="v29.0", category=FutureWarning) 

2072class WriteForcedSourceTableConfig(pipeBase.PipelineTaskConfig, 

2073 pipelineConnections=WriteForcedSourceTableConnections): 

2074 key = lsst.pex.config.Field( 

2075 doc="Column on which to join the two input tables on and make the primary key of the output", 

2076 dtype=str, 

2077 default="objectId", 

2078 ) 

2079 

2080 

2081@deprecated(reason="This task is replaced by lsst.pipe.tasks.ForcedPhotCcdTask. " 

2082 "This task will be removed after v30.", 

2083 version="v29.0", category=FutureWarning) 

2084class WriteForcedSourceTableTask(pipeBase.PipelineTask): 

2085 """Merge and convert per-detector forced source catalogs to DataFrame Parquet format. 

2086 

2087 Because the predecessor ForcedPhotCcdTask operates per-detector, 

2088 per-tract, (i.e., it has tract in its dimensions), detectors 

2089 on the tract boundary may have multiple forced source catalogs. 

2090 

2091 The successor task TransformForcedSourceTable runs per-patch 

2092 and temporally-aggregates overlapping mergedForcedSource catalogs from all 

2093 available multiple epochs. 

2094 """ 

2095 _DefaultName = "writeForcedSourceTable" 

2096 ConfigClass = WriteForcedSourceTableConfig 

2097 

2098 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

2099 inputs = butlerQC.get(inputRefs) 

2100 inputs["visit"] = butlerQC.quantum.dataId["visit"] 

2101 inputs["detector"] = butlerQC.quantum.dataId["detector"] 

2102 inputs["band"] = butlerQC.quantum.dataId["band"] 

2103 outputs = self.run(**inputs) 

2104 butlerQC.put(outputs, outputRefs) 

2105 

2106 def run(self, inputCatalog, inputCatalogDiff, visit, detector, band=None): 

2107 dfs = [] 

2108 for table, dataset, in zip((inputCatalog, inputCatalogDiff), ("calexp", "diff")): 

2109 df = table.asAstropy().to_pandas().set_index(self.config.key, drop=False) 

2110 df = df.reindex(sorted(df.columns), axis=1) 

2111 df["visit"] = visit 

2112 # int16 instead of uint8 because databases don't like unsigned bytes. 

2113 df["detector"] = np.int16(detector) 

2114 df["band"] = band if band else pd.NA 

2115 df.columns = pd.MultiIndex.from_tuples([(dataset, c) for c in df.columns], 

2116 names=("dataset", "column")) 

2117 

2118 dfs.append(df) 

2119 

2120 outputCatalog = functools.reduce(lambda d1, d2: d1.join(d2), dfs) 

2121 return pipeBase.Struct(outputCatalog=outputCatalog) 

2122 

2123 

2124class TransformForcedSourceTableConnections(pipeBase.PipelineTaskConnections, 

2125 dimensions=("instrument", "skymap", "patch", "tract")): 

2126 

2127 inputCatalogs = connectionTypes.Input( 

2128 doc="DataFrames of merged ForcedSources produced by WriteForcedSourceTableTask", 

2129 name="mergedForcedSource", 

2130 storageClass="DataFrame", 

2131 dimensions=("instrument", "visit", "detector", "skymap", "tract"), 

2132 multiple=True, 

2133 deferLoad=True 

2134 ) 

2135 referenceCatalog = connectionTypes.Input( 

2136 doc="Reference catalog which was used to seed the forcedPhot. Columns " 

2137 "objectId, detect_isPrimary, detect_isTractInner, detect_isPatchInner " 

2138 "are expected.", 

2139 name="objectTable", 

2140 storageClass="DataFrame", 

2141 dimensions=("tract", "patch", "skymap"), 

2142 deferLoad=True 

2143 ) 

2144 outputCatalog = connectionTypes.Output( 

2145 doc="Narrower, temporally-aggregated, per-patch ForcedSource Table transformed and converted per a " 

2146 "specified set of functors", 

2147 name="forcedSourceTable", 

2148 storageClass="ArrowAstropy", 

2149 dimensions=("tract", "patch", "skymap") 

2150 ) 

2151 

2152 

2153class TransformForcedSourceTableConfig(TransformCatalogBaseConfig, 

2154 pipelineConnections=TransformForcedSourceTableConnections): 

2155 referenceColumns = pexConfig.ListField( 

2156 dtype=str, 

2157 default=["detect_isPrimary", "detect_isTractInner", "detect_isPatchInner"], 

2158 optional=True, 

2159 doc="Columns to pull from reference catalog", 

2160 ) 

2161 keyRef = lsst.pex.config.Field( 

2162 doc="Column on which to join the two input tables on and make the primary key of the output", 

2163 dtype=str, 

2164 default="objectId", 

2165 ) 

2166 key = lsst.pex.config.Field( 

2167 doc="Rename the output DataFrame index to this name", 

2168 dtype=str, 

2169 default="forcedSourceId", 

2170 ) 

2171 

2172 def setDefaults(self): 

2173 super().setDefaults() 

2174 self.functorFile = os.path.join("$PIPE_TASKS_DIR", "schemas", "ForcedSource.yaml") 

2175 self.columnsFromDataId = ["tract", "patch"] 

2176 

2177 

2178class TransformForcedSourceTableTask(TransformCatalogBaseTask): 

2179 """Transform/standardize a ForcedSource catalog 

2180 

2181 Transforms each wide, per-detector forcedSource DataFrame per the 

2182 specification file (per-camera defaults found in ForcedSource.yaml). 

2183 All epochs that overlap the patch are aggregated into one per-patch 

2184 narrow-DataFrame file. 

2185 

2186 No de-duplication of rows is performed. Duplicate resolutions flags are 

2187 pulled in from the referenceCatalog: `detect_isPrimary`, 

2188 `detect_isTractInner`,`detect_isPatchInner`, so that user may de-duplicate 

2189 for analysis or compare duplicates for QA. 

2190 

2191 The resulting table includes multiple bands. Epochs (MJDs) and other useful 

2192 per-visit rows can be retreived by joining with the CcdVisitTable on 

2193 ccdVisitId. 

2194 """ 

2195 _DefaultName = "transformForcedSourceTable" 

2196 ConfigClass = TransformForcedSourceTableConfig 

2197 

2198 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

2199 inputs = butlerQC.get(inputRefs) 

2200 if self.funcs is None: 

2201 raise ValueError("config.functorFile is None. " 

2202 "Must be a valid path to yaml in order to run Task as a PipelineTask.") 

2203 outputs = self.run(inputs["inputCatalogs"], inputs["referenceCatalog"], funcs=self.funcs, 

2204 dataId=dict(outputRefs.outputCatalog.dataId.mapping)) 

2205 

2206 butlerQC.put(outputs, outputRefs) 

2207 

2208 def run(self, inputCatalogs, referenceCatalog, funcs=None, dataId=None, band=None): 

2209 dfs = [] 

2210 refColumns = list(self.config.referenceColumns) 

2211 refColumns.append(self.config.keyRef) 

2212 try: 

2213 ref = referenceCatalog.get(parameters={"columns": refColumns}) 

2214 except ValueError: 

2215 raise NoWorkFound(f"No objects to run forced photometry for {dataId}.") 

2216 if ref.index.name != self.config.keyRef: 

2217 # If the DataFrame we loaded was originally written as some other 

2218 # Parquet type, it probably doesn't have the index set. If it was 

2219 # written as a DataFrame, the index should already be set and 

2220 # trying to set it again would be an error, since it doens't exist 

2221 # as a regular column anymore. 

2222 ref.set_index(self.config.keyRef, inplace=True) 

2223 self.log.info("Aggregating %s input catalogs" % (len(inputCatalogs))) 

2224 for handle in inputCatalogs: 

2225 result = self.transform(None, handle, funcs, dataId) 

2226 # Filter for only rows that were detected on (overlap) the patch 

2227 dfs.append(result.df.join(ref, how="inner")) 

2228 

2229 outputCatalog = pd.concat(dfs) 

2230 

2231 if outputCatalog.empty: 

2232 raise NoWorkFound(f"No forced photometry rows for {dataId}.") 

2233 

2234 # Now that we are done joining on config.keyRef 

2235 # Change index to config.key by 

2236 outputCatalog.index.rename(self.config.keyRef, inplace=True) 

2237 # Add config.keyRef to the column list 

2238 outputCatalog.reset_index(inplace=True) 

2239 

2240 if "forcedSourceId" in outputCatalog.columns: 

2241 # Set the forcedSourceId to the index. This is specified in the 

2242 # ForcedSource.yaml 

2243 outputCatalog.set_index("forcedSourceId", inplace=True, verify_integrity=True) 

2244 # Rename it to the config.key 

2245 outputCatalog.index.rename(self.config.key, inplace=True) 

2246 

2247 self.log.info("Made a table of %d columns and %d rows", 

2248 len(outputCatalog.columns), len(outputCatalog)) 

2249 return pipeBase.Struct(outputCatalog=pandas_to_astropy(outputCatalog)) 

2250 

2251 

2252class ConsolidateTractConnections(pipeBase.PipelineTaskConnections, 

2253 defaultTemplates={"catalogType": ""}, 

2254 dimensions=("instrument", "tract")): 

2255 inputCatalogs = connectionTypes.Input( 

2256 doc="Input per-patch DataFrame Tables to be concatenated", 

2257 name="{catalogType}ForcedSourceTable", 

2258 storageClass="DataFrame", 

2259 dimensions=("tract", "patch", "skymap"), 

2260 multiple=True, 

2261 ) 

2262 

2263 outputCatalog = connectionTypes.Output( 

2264 doc="Output per-tract concatenation of DataFrame Tables", 

2265 name="{catalogType}ForcedSourceTable_tract", 

2266 storageClass="DataFrame", 

2267 dimensions=("tract", "skymap"), 

2268 ) 

2269 

2270 

2271class ConsolidateTractConfig(pipeBase.PipelineTaskConfig, 

2272 pipelineConnections=ConsolidateTractConnections): 

2273 

2274 doUseSchema = pexConfig.Field( 

2275 dtype=bool, 

2276 default=False, 

2277 doc="Use an existing schema to coerce the data types of the output columns." 

2278 ) 

2279 schemaDir = pexConfig.Field( 

2280 dtype=str, 

2281 doc="Path to the directory containing schema definitions.", 

2282 default=os.path.join("${SDM_SCHEMAS_DIR}", 

2283 "yml"), 

2284 optional=True, 

2285 ) 

2286 schemaFile = pexConfig.Field( 

2287 dtype=str, 

2288 doc="Yaml file specifying the schema of the output catalog.", 

2289 optional=True, 

2290 ) 

2291 tableName = pexConfig.Field( 

2292 dtype=str, 

2293 doc="Name of the table in the schema file to read.", 

2294 optional=True, 

2295 ) 

2296 

2297 

2298class ConsolidateTractTask(pipeBase.PipelineTask): 

2299 """Concatenate any per-patch, dataframe list into a single 

2300 per-tract DataFrame. 

2301 """ 

2302 _DefaultName = "ConsolidateTract" 

2303 ConfigClass = ConsolidateTractConfig 

2304 

2305 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

2306 inputs = butlerQC.get(inputRefs) 

2307 # Not checking at least one inputCatalog exists because that'd be an 

2308 # empty QG. 

2309 self.log.info("Concatenating %s per-patch %s Tables", 

2310 len(inputs["inputCatalogs"]), 

2311 inputRefs.inputCatalogs[0].datasetType.name) 

2312 df = pd.concat(inputs["inputCatalogs"]) 

2313 if self.config.doUseSchema: 

2314 schemaFile = os.path.join(self.config.schemaDir, self.config.schemaFile) 

2315 schema = readSdmSchemaFile(schemaFile) 

2316 df = convertDataFrameToSdmSchema(schema, df, tableName=self.config.tableName) 

2317 butlerQC.put(pipeBase.Struct(outputCatalog=df), outputRefs) 

2318 

2319 

2320class ConsolidateParentTractConnections( 

2321 pipeBase.PipelineTaskConnections, 

2322 dimensions=("instrument", "tract") 

2323): 

2324 inputCatalogs = connectionTypes.Input( 

2325 doc="Parents of the deblended objects", 

2326 name="object_parent_patch", 

2327 storageClass="SourceCatalog", 

2328 dimensions=("tract", "patch", "skymap"), 

2329 multiple=True, 

2330 ) 

2331 

2332 outputCatalog = connectionTypes.Output( 

2333 doc="Output per-tract concatenation of DataFrame Tables", 

2334 name="object_parent", 

2335 storageClass="ArrowAstropy", 

2336 dimensions=("tract", "skymap"), 

2337 ) 

2338 

2339 

2340class ConsolidateParentTractConfig( 

2341 pipeBase.PipelineTaskConfig, 

2342 pipelineConnections=ConsolidateParentTractConnections, 

2343): 

2344 pass 

2345 

2346 

2347class ConsolidateParentTractTask(pipeBase.PipelineTask): 

2348 """Concatenate any per-patch, dataframe list into a single 

2349 per-tract DataFrame. 

2350 """ 

2351 _DefaultName = "ConsolidateTract" 

2352 ConfigClass = ConsolidateParentTractConfig 

2353 

2354 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

2355 self.log.info("Concatenating %s per-patch %s Tables", 

2356 len(inputRefs.inputCatalogs), 

2357 inputRefs.inputCatalogs[0].datasetType.name) 

2358 

2359 tables = [] 

2360 for ref in inputRefs.inputCatalogs: 

2361 catalog = butlerQC.get(ref) 

2362 table = catalog.asAstropy() 

2363 

2364 # Rename columns 

2365 table.rename_column("id", "objectId") 

2366 table.rename_column("parent", "parentObjectId") 

2367 table.rename_column("merge_peak_sky", "sky_object") 

2368 

2369 # Add tract and patch columns 

2370 table["tract"] = ref.dataId["tract"] 

2371 table["patch"] = ref.dataId["patch"] 

2372 

2373 tables.append(table) 

2374 outputTable = astropy.table.vstack(tables, join_type="exact") 

2375 butlerQC.put(pipeBase.Struct(outputCatalog=outputTable), outputRefs)