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

901 statements  

« prev     ^ index     » next       coverage.py v7.16.2, created at 2026-10-24 09:42 +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 # This used to be 0.25 but was multiplied by 2*ln(2) to compensate 

1247 # for the bug fixed in DM-56611 

1248 default=0.34657359028, 

1249 ) 

1250 min_n_good_to_shift_flux_ratio = pexConfig.Field[int]( 

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

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

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

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

1255 default=0, 

1256 ) 

1257 

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

1259 size_column = self.size_column 

1260 size_model = np.sqrt( 

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

1262 ) 

1263 small = size_model < self.max_reff_compact 

1264 n_obj = len(table) 

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

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

1267 output = {} 

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

1269 if len(input_bands) > 1: 

1270 flux_psf, fluxerr_psf_sq, flux_model, fluxerr_model_sq = ( 

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

1272 for input_band in input_bands: 

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

1274 table, band=input_band) 

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

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

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

1278 flux_psf[good] += flux_psf_b[good] 

1279 flux_model[good] += flux_model_b[good] 

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

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

1282 fluxerr_psf = np.sqrt(fluxerr_psf_sq) 

1283 fluxerr_model = np.sqrt(fluxerr_model_sq) 

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

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

1286 else: 

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

1288 table, band=input_bands[0]) 

1289 

1290 psf_sn = flux_psf/fluxerr_psf 

1291 flux_ratio = np.array(flux_psf / flux_model) 

1292 

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

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

1295 # Attempt to correct any flux-independent systematic offset 

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

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

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

1299 

1300 flux_ratio_err = np.sqrt( 

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

1302 ) 

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

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

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

1306 # Make it sigmoid-like with a stretch 

1307 stretch = self.fluxerr_stretch 

1308 extendedness *= stretch 

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

1310 

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

1312 output[column_out] = extendedness 

1313 return output 

1314 

1315 

1316class ConsolidateObjectTableConnections(pipeBase.PipelineTaskConnections, 

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

1318 inputCatalogs = connectionTypes.Input( 

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

1320 name="objectTable", 

1321 storageClass="ArrowAstropy", 

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

1323 multiple=True, 

1324 deferLoad=True, 

1325 ) 

1326 outputCatalog = connectionTypes.Output( 

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

1328 name="objectTable_tract", 

1329 storageClass="ArrowAstropy", 

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

1331 ) 

1332 

1333 

1334class ConsolidateObjectTableConfig(pipeBase.PipelineTaskConfig, 

1335 pipelineConnections=ConsolidateObjectTableConnections): 

1336 actions = ConfigurableActionStructField[ComputeColumnsAction]( 

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

1338 ) 

1339 coaddName = pexConfig.Field[str]( 

1340 default="deep", 

1341 doc="Name of coadd" 

1342 ) 

1343 

1344 def setDefaults(self): 

1345 super().setDefaults() 

1346 self.actions.extendedness = ModelExtendednessColumnAction() 

1347 

1348 

1349class ConsolidateObjectTableTask(pipeBase.PipelineTask): 

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

1351 

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

1353 """ 

1354 _DefaultName = "consolidateObjectTable" 

1355 ConfigClass = ConsolidateObjectTableConfig 

1356 

1357 inputDataset = "objectTable" 

1358 outputDataset = "objectTable_tract" 

1359 

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

1361 inputs = butlerQC.get(inputRefs) 

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

1363 len(inputs["inputCatalogs"])) 

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

1365 for action in self.config.actions: 

1366 computed = action(table) 

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

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

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

1370 

1371 

1372class TransformSourceTableConnections(pipeBase.PipelineTaskConnections, 

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

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

1375 

1376 inputCatalog = connectionTypes.Input( 

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

1378 name="{catalogType}source", 

1379 storageClass="DataFrame", 

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

1381 deferLoad=True 

1382 ) 

1383 outputCatalog = connectionTypes.Output( 

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

1385 "specified set of functors", 

1386 name="{catalogType}sourceTable", 

1387 storageClass="ArrowAstropy", 

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

1389 ) 

1390 

1391 

1392class TransformSourceTableConfig(TransformCatalogBaseConfig, 

1393 pipelineConnections=TransformSourceTableConnections): 

1394 

1395 def setDefaults(self): 

1396 super().setDefaults() 

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

1398 self.primaryKey = "sourceId" 

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

1400 

1401 

1402class TransformSourceTableTask(TransformCatalogBaseTask): 

1403 """Transform/standardize a source catalog 

1404 """ 

1405 _DefaultName = "transformSourceTable" 

1406 ConfigClass = TransformSourceTableConfig 

1407 

1408 

1409class ConsolidateVisitSummaryConnections(pipeBase.PipelineTaskConnections, 

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

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

1412 camera = connectionTypes.PrerequisiteInput( 

1413 doc="Camera geometry.", 

1414 name="camera", 

1415 dimensions=("instrument",), 

1416 storageClass="Camera", 

1417 isCalibration=True, 

1418 ) 

1419 calexp = connectionTypes.Input( 

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

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

1422 name="calexp", 

1423 storageClass="ExposureF", 

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

1425 deferLoad=True, 

1426 multiple=True, 

1427 ) 

1428 visitSummary = connectionTypes.Output( 

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

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

1431 "detector."), 

1432 name="visitSummary", 

1433 storageClass="ExposureCatalog", 

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

1435 ) 

1436 visitSummarySchema = connectionTypes.InitOutput( 

1437 doc="Schema of the visitSummary catalog", 

1438 name="visitSummary_schema", 

1439 storageClass="ExposureCatalog", 

1440 ) 

1441 # RefitPointingTask options and connections use snake_case for consistency 

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

1443 # important than consistently using camelCase within this task. 

1444 visit_geometry = connectionTypes.Output( 

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

1446 name="visit_geometry", 

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

1448 storageClass="VisitGeometry", 

1449 ) 

1450 

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

1452 if self.config.do_refit_pointing: 

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

1454 else: 

1455 del self.camera 

1456 if not self.config.do_write_visit_geometry: 

1457 del self.visit_geometry 

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

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

1460 

1461 

1462class ConsolidateVisitSummaryConfig(pipeBase.PipelineTaskConfig, 

1463 pipelineConnections=ConsolidateVisitSummaryConnections): 

1464 """Config for ConsolidateVisitSummaryTask""" 

1465 

1466 full = pexConfig.Field( 

1467 "Whether to propate all exposure components. " 

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

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

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

1471 dtype=bool, 

1472 default=False, 

1473 ) 

1474 refitPointing = pexConfig.ConfigurableField( 

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

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

1477 target=RefitPointingTask, 

1478 ) 

1479 cameraDimensions = pexConfig.ListField( 

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

1481 dtype=str, 

1482 default=["instrument"], 

1483 ) 

1484 do_refit_pointing = pexConfig.Field( 

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

1486 dtype=bool, 

1487 default=True, 

1488 ) 

1489 do_write_visit_geometry = pexConfig.Field( 

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

1491 dtype=bool, 

1492 default=True, 

1493 ) 

1494 input_image_type = pexConfig.ChoiceField[str]( 

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

1496 allowed={ 

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

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

1499 }, 

1500 optional=False, 

1501 default="legacy", 

1502 ) 

1503 

1504 def validate(self): 

1505 super().validate() 

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

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

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

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

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

1511 

1512 def setDefaults(self): 

1513 super().setDefaults() 

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

1515 # UpdateVisitSummaryTask invocation of RefitPointingTask. 

1516 self.refitPointing.schema_prefix = "preliminary_" 

1517 

1518 

1519class ConsolidateVisitSummaryTask(pipeBase.PipelineTask): 

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

1521 

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

1523 single visit into an exposure catalog: 

1524 - The visitInfo. 

1525 - The wcs. 

1526 - The photoCalib. 

1527 - The physical_filter and band (if available). 

1528 - The PSF model. 

1529 - The aperture correction map. 

1530 - The transmission curve. 

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

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

1533 

1534 Tests for this task are performed in ci_hsc_gen3. 

1535 """ 

1536 _DefaultName = "consolidateVisitSummary" 

1537 ConfigClass = ConsolidateVisitSummaryConfig 

1538 

1539 def __init__(self, **kwargs): 

1540 super().__init__(**kwargs) 

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

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

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

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

1545 ExposureSummaryStats.update_schema(self.schema) 

1546 if self.config.do_refit_pointing: 

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

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

1549 

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

1551 handles = butlerQC.get(inputRefs.calexp) 

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

1553 

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

1555 len(handles), visit) 

1556 

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

1558 

1559 result = pipeBase.Struct() 

1560 try: 

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

1562 except pipeBase.AlgorithmError as e: 

1563 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

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

1565 ) 

1566 butlerQC.put(result, outputRefs) 

1567 raise error from e 

1568 

1569 butlerQC.put(result, outputRefs) 

1570 

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

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

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

1574 and other visit metadata. 

1575 

1576 Parameters 

1577 ---------- 

1578 visit : `int` 

1579 Visit identification number. 

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

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

1582 `lsst.images.VisitImage` datasets depending on 

1583 ``config.input_image_type``. 

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

1585 Camera geometry. Required if and only if 

1586 ``do_refit_pointing=True``. 

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

1588 Output struct to modify in-place. 

1589 

1590 Returns 

1591 ------- 

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

1593 Struct with the following attributes: 

1594 

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

1596 catalog with per-detector summary information. 

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

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

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

1600 """ 

1601 cat = afwTable.ExposureCatalog(self.schema) 

1602 cat.resize(len(handles)) 

1603 

1604 cat["visit"] = visit 

1605 

1606 # These are common for all detectors and should not 

1607 # be read in more than once. 

1608 visitInfo = None 

1609 filterLabel = None 

1610 

1611 for i, dataRef in enumerate(handles): 

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

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

1614 else: 

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

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

1617 visitInfo = components.visitInfo 

1618 filterLabel = components.filterLabel 

1619 summaryStats = components.summaryStats 

1620 

1621 rec = cat[i] 

1622 rec.setBBox(components.bbox) 

1623 rec.setVisitInfo(visitInfo) 

1624 rec.setWcs(components.wcs) 

1625 rec.setPhotoCalib(components.photoCalib) 

1626 rec.setValidPolygon(components.validPolygon) 

1627 

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

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

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

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

1632 

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

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

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

1636 summaryStats.update_record(rec) 

1637 

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

1639 raise pipeBase.NoWorkFound( 

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

1641 ) 

1642 

1643 metadata = dafBase.PropertyList() 

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

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

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

1647 cat.setMetadata(metadata) 

1648 

1649 cat.sort() 

1650 

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

1652 result = pipeBase.Struct() 

1653 result.visitSummary = cat 

1654 

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

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

1657 result.visit_geometry = refitPointingResult.regions 

1658 

1659 return result 

1660 

1661 @staticmethod 

1662 def _readLegacyComponents(handle, readVisitInfo): 

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

1664 `lsst.afw.image.Exposure`. 

1665 

1666 Parameters 

1667 ---------- 

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

1669 Handle to the exposure. 

1670 readVisitInfo : `bool` 

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

1672 for every detector in a visit. 

1673 

1674 Returns 

1675 ------- 

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

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

1678 ``readVisitInfo``), ``summaryStats``, ``wcs``, ``photoCalib``, 

1679 ``bbox`` and ``validPolygon`` attributes. 

1680 """ 

1681 return pipeBase.Struct( 

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

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

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

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

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

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

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

1689 ) 

1690 

1691 @staticmethod 

1692 def _readFutureComponents(handle, readVisitInfo): 

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

1694 and convert it to legacy types. 

1695 

1696 Parameters 

1697 ---------- 

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

1699 Handle to the visit image. 

1700 readVisitInfo : `bool` 

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

1702 for every detector in a visit. 

1703 

1704 Returns 

1705 ------- 

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

1707 Struct with the same attributes as `_readLegacyComponents`. 

1708 """ 

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

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

1711 names.append("obs_info") 

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

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

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

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

1716 else: 

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

1718 visitInfo = None 

1719 filterLabel = None 

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

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

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

1723 bounds = read["bounds"] 

1724 return pipeBase.Struct( 

1725 visitInfo=visitInfo, 

1726 filterLabel=filterLabel, 

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

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

1729 photoCalib=photoCalib, 

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

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

1732 ) 

1733 

1734 

1735class ConsolidateSourceTableConnections(pipeBase.PipelineTaskConnections, 

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

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

1738 inputCatalogs = connectionTypes.Input( 

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

1740 name="{catalogType}sourceTable", 

1741 storageClass="ArrowAstropy", 

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

1743 multiple=True, 

1744 deferLoad=True, 

1745 ) 

1746 outputCatalog = connectionTypes.Output( 

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

1748 name="{catalogType}sourceTable_visit", 

1749 storageClass="ArrowAstropy", 

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

1751 ) 

1752 

1753 

1754class ConsolidateSourceTableConfig(pipeBase.PipelineTaskConfig, 

1755 pipelineConnections=ConsolidateSourceTableConnections): 

1756 pass 

1757 

1758 

1759class ConsolidateSourceTableTask(pipeBase.PipelineTask): 

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

1761 """ 

1762 _DefaultName = "consolidateSourceTable" 

1763 ConfigClass = ConsolidateSourceTableConfig 

1764 

1765 inputDataset = "sourceTable" 

1766 outputDataset = "sourceTable_visit" 

1767 

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

1769 # Docstring inherited. 

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

1771 detectorOrder.sort() 

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

1773 inputs = butlerQC.get(inputRefs) 

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

1775 len(inputs["inputCatalogs"])) 

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

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

1778 

1779 

1780class MakeCcdVisitTableConnections(pipeBase.PipelineTaskConnections, 

1781 dimensions=("instrument",), 

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

1783 visitSummaryRefs = connectionTypes.Input( 

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

1785 name="finalVisitSummary", 

1786 storageClass="ExposureCatalog", 

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

1788 multiple=True, 

1789 deferLoad=True, 

1790 ) 

1791 outputCatalog = connectionTypes.Output( 

1792 doc="CCD and Visit metadata table", 

1793 name="ccdVisitTable", 

1794 storageClass="ArrowAstropy", 

1795 dimensions=("instrument",) 

1796 ) 

1797 

1798 

1799class MakeCcdVisitTableConfig(pipeBase.PipelineTaskConfig, 

1800 pipelineConnections=MakeCcdVisitTableConnections): 

1801 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

1802 

1803 

1804class MakeCcdVisitTableTask(pipeBase.PipelineTask): 

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

1806 """ 

1807 _DefaultName = "makeCcdVisitTable" 

1808 ConfigClass = MakeCcdVisitTableConfig 

1809 

1810 def run(self, visitSummaryRefs): 

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

1812 

1813 Parameters 

1814 ---------- 

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

1816 List of DeferredDatasetHandles pointing to exposure catalogs with 

1817 per-detector summary information. 

1818 

1819 Returns 

1820 ------- 

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

1822 Results struct with attribute: 

1823 

1824 ``outputCatalog`` 

1825 Catalog of ccd and visit information. 

1826 """ 

1827 ccdEntries = [] 

1828 for visitSummaryRef in visitSummaryRefs: 

1829 visitSummary = visitSummaryRef.get() 

1830 if not visitSummary: 

1831 continue 

1832 visitInfo = visitSummary[0].getVisitInfo() 

1833 

1834 # Strip provenance to prevent merge confusion. 

1835 strip_provenance_from_fits_header(visitSummary.metadata) 

1836 

1837 ccdEntry = {} 

1838 summaryTable = visitSummary.asAstropy() 

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

1840 "pixelScale", "zenithDistance", 

1841 "expTime", "zeroPoint", "psfSigma", 

1842 "psfArea", "psfIxx", "psfIyy", "psfIxy", 

1843 "skyBg", "skyNoise", 

1844 "psfAdaptiveThresholdValue", "psfAdaptiveIncludeThresholdMultiplier", 

1845 "nShapeletsStar", "shapeletsOnlyIqScore", "shapeletsIqScore", 

1846 "centroidDiffShapeletsVsSlotMedian", 

1847 "shapeletsStarEMedian", "shapeletsStarUnNormalizedEMedian", 

1848 "refCatSourceDensity", "astromOffsetMean", "astromOffsetStd", 

1849 "nPsfStar", "psfStarDeltaE1Median", "psfStarDeltaE2Median", 

1850 "psfStarDeltaE1Scatter", "psfStarDeltaE2Scatter", 

1851 "psfStarDeltaSizeMedian", "psfStarDeltaSizeScatter", 

1852 "psfStarScaledDeltaSizeScatter", "psfTraceRadiusDelta", 

1853 "psfApFluxDelta", "psfApCorrSigmaScaledDelta", 

1854 "maxDistToNearestPsf", "starEMedian", "starUnNormalizedEMedian", 

1855 "starComa1Median", "starComa2Median", 

1856 "starTrefoil1Median", "starTrefoil2Median", 

1857 "starKurtosisMedian", "starE41Median", "starE42Median", 

1858 "effTime", "effTimePsfSigmaScale", 

1859 "effTimeSkyBgScale", "effTimeZeroPointScale", 

1860 "magLim"] 

1861 ccdEntry = summaryTable[selectColumns] 

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

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

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

1865 # table. 

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

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

1868 

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

1870 # compatibility. To be removed after September 2023. 

1871 ccdEntry["decl"] = ccdEntry["dec"] 

1872 

1873 ccdEntry["ccdVisitId"] = [ 

1874 self.config.idGenerator.apply( 

1875 visitSummaryRef.dataId, 

1876 detector=detector_id, 

1877 is_exposure=False, 

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

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

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

1881 # use case for IdGenerator. This usage for a 

1882 # summary table is rare. 

1883 for detector_id in summaryTable["id"] 

1884 ] 

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

1886 ccdEntry["seeing"] = ( 

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

1888 ) 

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

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

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

1892 expTime = visitInfo.getExposureTime() 

1893 ccdEntry["obsStart"] = ( 

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

1895 ) 

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

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

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

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

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

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

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

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

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

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

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

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

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

1909 # These columns are only added on certain updateVisitSummary 

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

1911 # convert to arcseconds here. 

1912 for inName, outName in { 

1913 "wcs_corner_max_offset": "wcsCornerMaxOffset", 

1914 "wcs_detector_pointing_residual": "wcsDetectorPointingResidual", 

1915 "wcs_visit_pointing_residual": "wcsVisitPointingResidual", 

1916 "preliminary_wcs_detector_pointing_residual": "wcsPreliminaryDetectorPointingResidual", 

1917 "preliminary_wcs_visit_pointing_residual": "wcsPreliminaryVisitPointingResidual", 

1918 "wcs_sip_delta_sky": "wcsSIPDeltaSky", 

1919 }.items(): 

1920 if inName in summaryTable.columns: 

1921 inCol = summaryTable[inName] 

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

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

1924 unit=astropy.units.arcsec 

1925 ) 

1926 if "wcs_sip_delta_pixel" in summaryTable.columns: 

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

1928 

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

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

1931 # values are actually wanted. 

1932 ccdEntries.append(ccdEntry) 

1933 

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

1935 return pipeBase.Struct(outputCatalog=outputCatalog) 

1936 

1937 

1938class MakeVisitTableConnections(pipeBase.PipelineTaskConnections, 

1939 dimensions=("instrument",), 

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

1941 visitSummaries = connectionTypes.Input( 

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

1943 name="finalVisitSummary", 

1944 storageClass="ExposureCatalog", 

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

1946 multiple=True, 

1947 deferLoad=True, 

1948 ) 

1949 outputCatalog = connectionTypes.Output( 

1950 doc="Visit metadata table", 

1951 name="visitTable", 

1952 storageClass="ArrowAstropy", 

1953 dimensions=("instrument",) 

1954 ) 

1955 

1956 

1957class MakeVisitTableConfig(pipeBase.PipelineTaskConfig, 

1958 pipelineConnections=MakeVisitTableConnections): 

1959 pass 

1960 

1961 

1962class MakeVisitTableTask(pipeBase.PipelineTask): 

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

1964 """ 

1965 _DefaultName = "makeVisitTable" 

1966 ConfigClass = MakeVisitTableConfig 

1967 

1968 def run(self, visitSummaries): 

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

1970 

1971 Parameters 

1972 ---------- 

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

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

1975 Returns 

1976 ------- 

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

1978 Results struct with attribute: 

1979 

1980 ``outputCatalog`` 

1981 Catalog of visit information. 

1982 """ 

1983 visitEntries = [] 

1984 for visitSummary in visitSummaries: 

1985 visitSummary = visitSummary.get() 

1986 if not visitSummary: 

1987 continue 

1988 visitRow = visitSummary[0] 

1989 visitInfo = visitRow.getVisitInfo() 

1990 

1991 visitEntry = {} 

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

1993 visitEntry["visit"] = visitRow["visit"] 

1994 visitEntry["physical_filter"] = visitRow["physical_filter"] 

1995 visitEntry["band"] = visitRow["band"] 

1996 raDec = visitInfo.getBoresightRaDec() 

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

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

1999 

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

2001 # compatibility. To be removed after September 2023. 

2002 visitEntry["decl"] = visitEntry["dec"] 

2003 

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

2005 azAlt = visitInfo.getBoresightAzAlt() 

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

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

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

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

2010 expTime = visitInfo.getExposureTime() 

2011 visitEntry["expTime"] = expTime 

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

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

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

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

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

2017 

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

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

2020 visitEntry["psfTE1e1"] = visitRow["psfTE1e1"] 

2021 visitEntry["psfTE1e2"] = visitRow["psfTE1e2"] 

2022 visitEntry["psfTE1ex"] = visitRow["psfTE1ex"] 

2023 visitEntry["psfTE2e1"] = visitRow["psfTE2e1"] 

2024 visitEntry["psfTE2e2"] = visitRow["psfTE2e2"] 

2025 visitEntry["psfTE2ex"] = visitRow["psfTE2ex"] 

2026 visitEntry["psfTE3e1"] = visitRow["psfTE3e1"] 

2027 visitEntry["psfTE3e2"] = visitRow["psfTE3e2"] 

2028 visitEntry["psfTE3ex"] = visitRow["psfTE3ex"] 

2029 visitEntry["psfTE4e1"] = visitRow["psfTE4e1"] 

2030 visitEntry["psfTE4e2"] = visitRow["psfTE4e2"] 

2031 visitEntry["psfTE4ex"] = visitRow["psfTE4ex"] 

2032 

2033 visitEntries.append(visitEntry) 

2034 

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

2036 # mirror1Temp, mirror2Temp, mirror3Temp, domeTemp, externalTemp, 

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

2038 

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

2040 return pipeBase.Struct(outputCatalog=outputCatalog) 

2041 

2042 

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

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

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

2046class WriteForcedSourceTableConnections(pipeBase.PipelineTaskConnections, 

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

2048 

2049 inputCatalog = connectionTypes.Input( 

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

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

2052 name="forced_src", 

2053 storageClass="SourceCatalog", 

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

2055 ) 

2056 inputCatalogDiff = connectionTypes.Input( 

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

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

2059 name="forced_diff", 

2060 storageClass="SourceCatalog", 

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

2062 ) 

2063 outputCatalog = connectionTypes.Output( 

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

2065 name="mergedForcedSource", 

2066 storageClass="DataFrame", 

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

2068 ) 

2069 

2070 

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

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

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

2074class WriteForcedSourceTableConfig(pipeBase.PipelineTaskConfig, 

2075 pipelineConnections=WriteForcedSourceTableConnections): 

2076 key = lsst.pex.config.Field( 

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

2078 dtype=str, 

2079 default="objectId", 

2080 ) 

2081 

2082 

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

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

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

2086class WriteForcedSourceTableTask(pipeBase.PipelineTask): 

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

2088 

2089 Because the predecessor ForcedPhotCcdTask operates per-detector, 

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

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

2092 

2093 The successor task TransformForcedSourceTable runs per-patch 

2094 and temporally-aggregates overlapping mergedForcedSource catalogs from all 

2095 available multiple epochs. 

2096 """ 

2097 _DefaultName = "writeForcedSourceTable" 

2098 ConfigClass = WriteForcedSourceTableConfig 

2099 

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

2101 inputs = butlerQC.get(inputRefs) 

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

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

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

2105 outputs = self.run(**inputs) 

2106 butlerQC.put(outputs, outputRefs) 

2107 

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

2109 dfs = [] 

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

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

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

2113 df["visit"] = visit 

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

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

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

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

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

2119 

2120 dfs.append(df) 

2121 

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

2123 return pipeBase.Struct(outputCatalog=outputCatalog) 

2124 

2125 

2126class TransformForcedSourceTableConnections(pipeBase.PipelineTaskConnections, 

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

2128 

2129 inputCatalogs = connectionTypes.Input( 

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

2131 name="mergedForcedSource", 

2132 storageClass="DataFrame", 

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

2134 multiple=True, 

2135 deferLoad=True 

2136 ) 

2137 referenceCatalog = connectionTypes.Input( 

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

2139 "objectId, detect_isPrimary, detect_isTractInner, detect_isPatchInner " 

2140 "are expected.", 

2141 name="objectTable", 

2142 storageClass="DataFrame", 

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

2144 deferLoad=True 

2145 ) 

2146 outputCatalog = connectionTypes.Output( 

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

2148 "specified set of functors", 

2149 name="forcedSourceTable", 

2150 storageClass="ArrowAstropy", 

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

2152 ) 

2153 

2154 

2155class TransformForcedSourceTableConfig(TransformCatalogBaseConfig, 

2156 pipelineConnections=TransformForcedSourceTableConnections): 

2157 referenceColumns = pexConfig.ListField( 

2158 dtype=str, 

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

2160 optional=True, 

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

2162 ) 

2163 keyRef = lsst.pex.config.Field( 

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

2165 dtype=str, 

2166 default="objectId", 

2167 ) 

2168 key = lsst.pex.config.Field( 

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

2170 dtype=str, 

2171 default="forcedSourceId", 

2172 ) 

2173 

2174 def setDefaults(self): 

2175 super().setDefaults() 

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

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

2178 

2179 

2180class TransformForcedSourceTableTask(TransformCatalogBaseTask): 

2181 """Transform/standardize a ForcedSource catalog 

2182 

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

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

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

2186 narrow-DataFrame file. 

2187 

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

2189 pulled in from the referenceCatalog: `detect_isPrimary`, 

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

2191 for analysis or compare duplicates for QA. 

2192 

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

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

2195 ccdVisitId. 

2196 """ 

2197 _DefaultName = "transformForcedSourceTable" 

2198 ConfigClass = TransformForcedSourceTableConfig 

2199 

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

2201 inputs = butlerQC.get(inputRefs) 

2202 if self.funcs is None: 

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

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

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

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

2207 

2208 butlerQC.put(outputs, outputRefs) 

2209 

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

2211 dfs = [] 

2212 refColumns = list(self.config.referenceColumns) 

2213 refColumns.append(self.config.keyRef) 

2214 try: 

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

2216 except ValueError: 

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

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

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

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

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

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

2223 # as a regular column anymore. 

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

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

2226 for handle in inputCatalogs: 

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

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

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

2230 

2231 outputCatalog = pd.concat(dfs) 

2232 

2233 if outputCatalog.empty: 

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

2235 

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

2237 # Change index to config.key by 

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

2239 # Add config.keyRef to the column list 

2240 outputCatalog.reset_index(inplace=True) 

2241 

2242 if "forcedSourceId" in outputCatalog.columns: 

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

2244 # ForcedSource.yaml 

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

2246 # Rename it to the config.key 

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

2248 

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

2250 len(outputCatalog.columns), len(outputCatalog)) 

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

2252 

2253 

2254class ConsolidateTractConnections(pipeBase.PipelineTaskConnections, 

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

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

2257 inputCatalogs = connectionTypes.Input( 

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

2259 name="{catalogType}ForcedSourceTable", 

2260 storageClass="DataFrame", 

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

2262 multiple=True, 

2263 ) 

2264 

2265 outputCatalog = connectionTypes.Output( 

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

2267 name="{catalogType}ForcedSourceTable_tract", 

2268 storageClass="DataFrame", 

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

2270 ) 

2271 

2272 

2273class ConsolidateTractConfig(pipeBase.PipelineTaskConfig, 

2274 pipelineConnections=ConsolidateTractConnections): 

2275 

2276 doUseSchema = pexConfig.Field( 

2277 dtype=bool, 

2278 default=False, 

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

2280 ) 

2281 schemaDir = pexConfig.Field( 

2282 dtype=str, 

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

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

2285 "yml"), 

2286 optional=True, 

2287 ) 

2288 schemaFile = pexConfig.Field( 

2289 dtype=str, 

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

2291 optional=True, 

2292 ) 

2293 tableName = pexConfig.Field( 

2294 dtype=str, 

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

2296 optional=True, 

2297 ) 

2298 

2299 

2300class ConsolidateTractTask(pipeBase.PipelineTask): 

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

2302 per-tract DataFrame. 

2303 """ 

2304 _DefaultName = "ConsolidateTract" 

2305 ConfigClass = ConsolidateTractConfig 

2306 

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

2308 inputs = butlerQC.get(inputRefs) 

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

2310 # empty QG. 

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

2312 len(inputs["inputCatalogs"]), 

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

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

2315 if self.config.doUseSchema: 

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

2317 schema = readSdmSchemaFile(schemaFile) 

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

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

2320 

2321 

2322class ConsolidateParentTractConnections( 

2323 pipeBase.PipelineTaskConnections, 

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

2325): 

2326 inputCatalogs = connectionTypes.Input( 

2327 doc="Parents of the deblended objects", 

2328 name="object_parent_patch", 

2329 storageClass="SourceCatalog", 

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

2331 multiple=True, 

2332 ) 

2333 

2334 outputCatalog = connectionTypes.Output( 

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

2336 name="object_parent", 

2337 storageClass="ArrowAstropy", 

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

2339 ) 

2340 

2341 

2342class ConsolidateParentTractConfig( 

2343 pipeBase.PipelineTaskConfig, 

2344 pipelineConnections=ConsolidateParentTractConnections, 

2345): 

2346 pass 

2347 

2348 

2349class ConsolidateParentTractTask(pipeBase.PipelineTask): 

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

2351 per-tract DataFrame. 

2352 """ 

2353 _DefaultName = "ConsolidateTract" 

2354 ConfigClass = ConsolidateParentTractConfig 

2355 

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

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

2358 len(inputRefs.inputCatalogs), 

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

2360 

2361 tables = [] 

2362 for ref in inputRefs.inputCatalogs: 

2363 catalog = butlerQC.get(ref) 

2364 table = catalog.asAstropy() 

2365 

2366 # Rename columns 

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

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

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

2370 

2371 # Add tract and patch columns 

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

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

2374 

2375 tables.append(table) 

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

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