Coverage for python / lsst / drp / tasks / metadetection_shear.py: 20%

1# This file is part of drp_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 

22from __future__ import annotations 

23 

24__all__ = ( 

25 "MetadetectionProcessingError", 

26 "MetadetectionShearConfig", 

27 "MetadetectionShearTask", 

28) 

29 

30from collections.abc import Collection, Mapping, Sequence 

31from itertools import product 

32from typing import Any, ClassVar 

33 

34import esutil as eu 

35import numpy as np 

36import pyarrow as pa 

37from metadetect.lsst.masking import apply_apodized_bright_masks_mbexp, apply_apodized_edge_masks_mbexp 

38from metadetect.lsst.metacal_exposures import STEP as SHEAR_STEP 

39from metadetect.lsst.metadetect import MetadetectTask 

40from metadetect.lsst.util import extract_multiband_coadd_data 

41 

42import lsst.pipe.base.connectionTypes as cT 

43from lsst.afw.image import ExposureF 

44from lsst.afw.table import SimpleCatalog 

45from lsst.cell_coadds import MultipleCellCoadd, StitchedCoadd 

46from lsst.daf.butler import DataCoordinate, DatasetRef 

47from lsst.meas.algorithms import LoadReferenceObjectsConfig, ReferenceObjectLoader 

48from lsst.meas.base import FullIdGenerator, SkyMapIdGeneratorConfig 

49from lsst.pex.config import ConfigField, ConfigurableField, Field, FieldValidationError, ListField 

50from lsst.pipe.base import ( 

51 AlgorithmError, 

52 InputQuantizedConnection, 

53 NoWorkFound, 

54 OutputQuantizedConnection, 

55 PipelineTask, 

56 PipelineTaskConfig, 

57 PipelineTaskConnections, 

58 QuantumContext, 

59 Struct, 

60) 

61from lsst.pipe.base.connectionTypes import BaseInput, Output 

62from lsst.skymap import Index2D 

63 

64 

65class MetadetectionProcessingError(AlgorithmError): 

66 """Exception raised when metadetection processing fails.""" 

67 

68 @property 

69 def metadata(self) -> dict: 

70 return {} 

71 

72 

73class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch"}): 

74 """Definitions of inputs and outputs for MetadetectionShearTask.""" 

75 

76 input_coadds = cT.Input( 

77 "deep_coadd_cell_predetection", 

78 storageClass="MultipleCellCoadd", 

79 doc="Per-band deep coadds.", 

80 multiple=True, 

81 dimensions={"patch", "band"}, 

82 ) 

83 

84 ref_cat = cT.PrerequisiteInput( 

85 doc="Reference catalog used to mask bright objects.", 

86 name="the_monster_20250219", 

87 storageClass="SimpleCatalog", 

88 dimensions=("skypix",), 

89 deferLoad=True, 

90 multiple=True, 

91 ) 

92 

93 metadetect_catalog = cT.Output( 

94 "object_shear_patch", 

95 storageClass="ArrowTable", 

96 doc="Output catalog with all quantities measured inside the metacalibration loop.", 

97 multiple=False, 

98 dimensions={"patch"}, 

99 ) 

100 

101 metadetect_schema = cT.InitOutput( 

102 "object_shear_schema", 

103 storageClass="ArrowSchema", 

104 doc="Schema of the output catalog.", 

105 ) 

106 

107 config: MetadetectionShearConfig 

108 

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

110 super().__init__(config=config) 

111 

112 if not config: 

113 return None 

114 

115 if not config.do_mask_bright_objects: 

116 del self.ref_cat 

117 

118 def adjustQuantum( 

119 self, 

120 inputs: dict[str, tuple[BaseInput, Collection[DatasetRef]]], 

121 outputs: dict[str, tuple[Output, Collection[DatasetRef]]], 

122 label: str, 

123 data_id: DataCoordinate, 

124 ) -> tuple[ 

125 Mapping[str, tuple[BaseInput, Collection[DatasetRef]]], 

126 Mapping[str, tuple[Output, Collection[DatasetRef]]], 

127 ]: 

128 # Docstring inherited. 

129 # This is a hook for customizing what is input and output to each 

130 # invocation of the task as early as possible, which we override here 

131 # to make sure we have exactly the required bands, no more, no less. 

132 connection, original_input_coadds = inputs["input_coadds"] 

133 bands_missing = set(self.config.photometry_bands) 

134 adjusted_input_coadds = [] 

135 for ref in original_input_coadds: 

136 if ref.dataId["band"] in self.config.photometry_bands: 

137 adjusted_input_coadds.append(ref) 

138 bands_missing.remove(ref.dataId["band"]) 

139 if missing_shear_bands := bands_missing.intersection(self.config.metadetect.shear_bands): 

140 raise NoWorkFound(f"Required bands {missing_shear_bands} not present for {label}@{data_id}).") 

141 adjusted_inputs = {"input_coadds": (connection, adjusted_input_coadds)} 

142 inputs.update(adjusted_inputs) 

143 super().adjustQuantum(inputs, outputs, label, data_id) 

144 return adjusted_inputs, {} 

145 

146 

147class MetadetectionShearConfig(PipelineTaskConfig, pipelineConnections=MetadetectionShearConnections): 

148 """Configuration definition for MetadetectionShearTask.""" 

149 

150 metadetect = ConfigurableField( 

151 target=MetadetectTask, 

152 doc="Configuration for metadetection.", 

153 ) 

154 

155 photometry_bands = ListField[str]( 

156 "Bands expected to be present. Cells with one or more of these bands " 

157 "missing will be skipped. Bands other than those listed here will " 

158 "not be processed.", 

159 default=["g", "r", "i", "z"], 

160 ) 

161 

162 do_mask_bright_objects = Field[bool]( 

163 doc="Mask bright objects in coadds?", 

164 default=False, 

165 ) 

166 

167 ref_loader = ConfigField( 

168 dtype=LoadReferenceObjectsConfig, 

169 doc="Reference object loader used for bright-object masking.", 

170 ) 

171 

172 ref_loader_filter_name = Field[str]( 

173 "Filter name from ref_loader used for bright-object masking.", 

174 default="monster_DES_r", 

175 ) 

176 

177 border = Field[int]( 

178 "Border to apply to single cell images", 

179 default=50, 

180 ) 

181 

182 id_generator = SkyMapIdGeneratorConfig.make_field() 

183 

184 def setDefaults(self): 

185 super().setDefaults() 

186 self.metadetect.shear_bands = ["r", "i", "z"] 

187 self.metadetect.metacal.types = ["noshear", "1p", "1m", "2p", "2m"] 

188 

189 def validate(self): 

190 super().validate() 

191 if (shear_bands := self.metadetect.shear_bands) is not None and not set(shear_bands).issubset( 

192 self.photometry_bands 

193 ): 

194 raise FieldValidationError( 

195 self.__class__.metadetect, 

196 self, 

197 "photometry_bands must be a list of bands that is a superset of metadetect.shear_bands", 

198 ) 

199 

200 

201class MetadetectionShearTask(PipelineTask): 

202 """A PipelineTask that measures shear using metadetection.""" 

203 

204 _DefaultName: ClassVar[str] = "metadetectionShear" 

205 ConfigClass: ClassVar[type[MetadetectionShearConfig]] = MetadetectionShearConfig 

206 

207 config: MetadetectionShearConfig 

208 

209 def __init__(self, *, initInputs: dict[str, Any] | None = None, **kwargs: Any): 

210 super().__init__(initInputs=initInputs, **kwargs) 

211 self.metadetect_schema = self.make_metadetect_schema(self.config) 

212 self.makeSubtask("metadetect") 

213 

214 @classmethod 

215 def make_metadetect_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: 

216 """Construct a PyArrow Schema for this task's main output catalog. 

217 

218 Parameters 

219 ---------- 

220 config : `MetadetectionShearConfig` 

221 Configuration that may be used to control details of the schema. 

222 

223 Returns 

224 ------- 

225 object_schema : `pyarrow.Schema` 

226 Schema for the object catalog produced by this task. Each field's 

227 metadata should include both a 'doc' entry and a 'unit' entry. 

228 """ 

229 pa_schema = pa.schema( 

230 [ 

231 # Fields from pipeline bookkeeping. 

232 pa.field( 

233 "shearObjectId", 

234 pa.uint64(), 

235 nullable=False, 

236 metadata={ 

237 "doc": ( 

238 "Unique identifier for a ShearObject, specific " 

239 "to a single metacalibration counterfactual image." 

240 ), 

241 "unit": "", 

242 }, 

243 ), 

244 pa.field( 

245 "tract", 

246 pa.uint64(), 

247 nullable=False, 

248 metadata={ 

249 "doc": "ID of the tract on which this measurement was made.", 

250 "unit": "", 

251 }, 

252 ), 

253 pa.field( 

254 "patch", 

255 pa.uint32(), 

256 nullable=False, 

257 metadata={ 

258 "doc": "ID of the patch within the tract on which this measurement was made.", 

259 "unit": "", 

260 }, 

261 ), 

262 pa.field( 

263 "cell_x", 

264 pa.uint32(), 

265 nullable=False, 

266 metadata={ 

267 "doc": "Column of the cell within the patch on which this measurement was made.", 

268 "unit": "", 

269 }, 

270 ), 

271 pa.field( 

272 "cell_y", 

273 pa.uint32(), 

274 nullable=False, 

275 metadata={ 

276 "doc": "Row of the cell within the patch on which this measurement was made.", 

277 "unit": "", 

278 }, 

279 ), 

280 # Fields from metadetection (generic). 

281 pa.field( 

282 "metaStep", 

283 pa.string(), 

284 nullable=False, 

285 metadata={ 

286 "doc": ( 

287 "Type of artificial shear applied to image. " 

288 "One of: 'ns', '1p', '1m', '2p', '2m'." 

289 ), 

290 "unit": "", 

291 }, 

292 ), 

293 pa.field( 

294 "stamp_flags", 

295 pa.uint32(), 

296 nullable=False, 

297 metadata={ 

298 "doc": "Flags for the stamp on which this measurement was made.", 

299 "unit": "", 

300 }, 

301 ), 

302 pa.field( 

303 "x", 

304 pa.float32(), 

305 nullable=False, 

306 metadata={ 

307 "doc": "Centroid (tract, x-axis) of the detected ShearObject.", 

308 "unit": "", 

309 }, 

310 ), 

311 pa.field( 

312 "y", 

313 pa.float32(), 

314 nullable=False, 

315 metadata={ 

316 "doc": "Centroid (tract, y-axis) of the detected ShearObject.", 

317 "unit": "", 

318 }, 

319 ), 

320 pa.field( 

321 "coord_ra", 

322 pa.float64(), 

323 nullable=False, 

324 metadata={ 

325 "doc": "Detected Right Ascension of the ShearObject.", 

326 "unit": "degrees", 

327 }, 

328 ), 

329 pa.field( 

330 "coord_dec", 

331 pa.float64(), 

332 nullable=False, 

333 metadata={ 

334 "doc": "Detected Declination of the ShearObject.", 

335 "unit": "degrees", 

336 }, 

337 ), 

338 # Original PSF measurements 

339 pa.field( 

340 "psfOriginal_flags", 

341 pa.uint32(), 

342 nullable=False, 

343 metadata={ 

344 "doc": "Flags for the original PSF measurement.", 

345 "unit": "", 

346 }, 

347 ), 

348 pa.field( 

349 "psfOriginal_e1", 

350 pa.float32(), 

351 nullable=False, 

352 metadata={ 

353 "doc": "Distortion-style e1 of the original PSF from adaptive moments.", 

354 "unit": "", 

355 }, 

356 ), 

357 pa.field( 

358 "psfOriginal_e2", 

359 pa.float32(), 

360 nullable=False, 

361 metadata={ 

362 "doc": "Distortion-style e2 of the original PSF from adaptive moments.", 

363 "unit": "", 

364 }, 

365 ), 

366 pa.field( 

367 "psfOriginal_T", 

368 pa.float32(), 

369 nullable=False, 

370 metadata={ 

371 "doc": "Trace (<x^2> + <y^2>) measurement of the original PSF from adaptive moments.", 

372 "unit": "arcseconds squared", 

373 }, 

374 ), 

375 pa.field( 

376 "bmask", 

377 pa.uint32(), 

378 nullable=False, 

379 metadata={ 

380 "doc": "`bmask` flags for the ShearObject", 

381 "unit": "", 

382 }, 

383 ), 

384 pa.field( 

385 "ormask", 

386 pa.int32(), 

387 nullable=False, 

388 metadata={ 

389 "doc": "`ored` mask flags for the ShearObject", 

390 "unit": "", 

391 }, 

392 ), 

393 pa.field( 

394 "mfrac", 

395 pa.float32(), 

396 nullable=False, 

397 metadata={ 

398 "doc": "Gaussian-weighted masked fraction for the ShearObject.", 

399 "unit": "", 

400 }, 

401 ), 

402 # Fields that come only from gauss algorithm. 

403 # Reconvolved PSF measurements (gauss) 

404 pa.field( 

405 "gauss_psfReconvolved_flags", 

406 pa.uint32(), 

407 nullable=False, 

408 metadata={ 

409 "doc": "Flags for reconvolved PSF (measured with gauss algorithm).", 

410 "unit": "", 

411 }, 

412 ), 

413 pa.field( 

414 "gauss_psfReconvolved_g1", 

415 pa.float32(), 

416 nullable=False, 

417 metadata={ 

418 "doc": "Reduced-shear g1 of the reconvolved PSF (measured with gauss algorithm).", 

419 "unit": "", 

420 }, 

421 ), 

422 pa.field( 

423 "gauss_psfReconvolved_g2", 

424 pa.float32(), 

425 nullable=False, 

426 metadata={ 

427 "doc": "Reduced-shear g2 of the reconvolved PSF (measured with gauss algorithm).", 

428 "unit": "", 

429 }, 

430 ), 

431 pa.field( 

432 "gauss_psfReconvolved_T", 

433 pa.float32(), 

434 nullable=False, 

435 metadata={ 

436 "doc": ( 

437 "Trace (<x^2> + <y^2>) of the reconvolved PSF (measured with gauss algorithm)." 

438 ), 

439 "unit": "arcseconds squared", 

440 }, 

441 ), 

442 # Object measurements (gauss algorithm). 

443 pa.field( 

444 "gauss_g1", 

445 pa.float32(), 

446 nullable=False, 

447 metadata={ 

448 "doc": ( 

449 "Reduced-shear g1 measurement of the ShearObject " 

450 "(measured with gauss algorithm)." 

451 ), 

452 "unit": "", 

453 }, 

454 ), 

455 pa.field( 

456 "gauss_g2", 

457 pa.float32(), 

458 nullable=False, 

459 metadata={ 

460 "doc": ( 

461 "Reduced-shear g2 measurement of the ShearObject " 

462 "(measured with gauss algorithm)." 

463 ), 

464 "unit": "", 

465 }, 

466 ), 

467 pa.field( 

468 "gauss_g1_g1_Cov", 

469 pa.float32(), 

470 nullable=False, 

471 metadata={ 

472 "doc": ( 

473 "Auto-covariance of g1 measurement of the ShearObject " 

474 "(measured with gauss algorithm)." 

475 ), 

476 "unit": "", 

477 }, 

478 ), 

479 pa.field( 

480 "gauss_g1_g2_Cov", 

481 pa.float32(), 

482 nullable=False, 

483 metadata={ 

484 "doc": ( 

485 "Cross-covariance of g1 and g2 measurement of the ShearObject " 

486 "(measured with gauss algorithm)." 

487 ), 

488 "unit": "", 

489 }, 

490 ), 

491 pa.field( 

492 "gauss_g2_g2_Cov", 

493 pa.float32(), 

494 nullable=False, 

495 metadata={ 

496 "doc": ( 

497 "Auto-covariance of g2 measurement of the ShearObject " 

498 "(measured with gauss algorithm)." 

499 ), 

500 "unit": "", 

501 }, 

502 ), 

503 ], 

504 metadata={ 

505 "shear_step": str(SHEAR_STEP), 

506 "shear_bands": "".join(sorted(config.metadetect.shear_bands)), 

507 }, 

508 ) 

509 

510 for alg_name in ("gauss", "pgauss"): 

511 pa_schema = pa_schema.append( 

512 pa.field( 

513 f"{alg_name}_snr", 

514 pa.float32(), 

515 nullable=False, 

516 metadata={ 

517 "doc": ( 

518 "Signal-to-noise ratio measure of the ShearObject " 

519 f"(measured with {alg_name} algorithm)." 

520 ), 

521 "unit": "", 

522 }, 

523 ), 

524 ) 

525 pa_schema = pa_schema.append( 

526 pa.field( 

527 f"{alg_name}_T", 

528 pa.float32(), 

529 nullable=False, 

530 metadata={ 

531 "doc": ( 

532 "Trace (<x^2> + <y^2>) measurement of the ShearObject " 

533 f"(measured with {alg_name} algorithm)." 

534 ), 

535 "unit": "arcseconds squared", 

536 }, 

537 ), 

538 ) 

539 pa_schema = pa_schema.append( 

540 pa.field( 

541 f"{alg_name}_TErr", 

542 pa.float32(), 

543 nullable=False, 

544 metadata={ 

545 "doc": ( 

546 "Uncertainty in the trace measurement of the ShearObject " 

547 f"(measured with {alg_name} algorithm)." 

548 ), 

549 "unit": "arcseconds squared", 

550 }, 

551 ), 

552 ) 

553 pa_schema = pa_schema.append( 

554 pa.field( 

555 f"{alg_name}_T_flags", 

556 pa.uint32(), 

557 nullable=False, 

558 metadata={ 

559 "doc": ( 

560 "Flags for the trace (<x^2> + <y^2>) measurement of the ShearObject " 

561 f"(measured with {alg_name} algorithm)." 

562 ), 

563 "unit": "", 

564 }, 

565 ), 

566 ) 

567 pa_schema = pa_schema.append( 

568 pa.field( 

569 f"{alg_name}_object_flags", 

570 pa.uint32(), 

571 nullable=False, 

572 metadata={ 

573 "doc": f"Flags for the ShearObject measurement (measured with {alg_name} algorithm).", 

574 "unit": "", 

575 }, 

576 ), 

577 ) 

578 pa_schema = pa_schema.append( 

579 pa.field( 

580 f"{alg_name}_flags", 

581 pa.uint32(), 

582 nullable=False, 

583 metadata={ 

584 "doc": f"Overall flags for {alg_name} measurement algorithm.", 

585 "unit": "", 

586 }, 

587 ), 

588 ) 

589 

590 # Per-band quantities, typically fluxes and associated quantites. 

591 for b in config.photometry_bands: 

592 pa_schema = pa_schema.append( 

593 pa.field( 

594 f"{b}_{alg_name}Flux_flags", 

595 pa.uint32(), 

596 nullable=False, 

597 metadata={ 

598 "doc": f"Flags set for flux in {b} band measured with {alg_name} algorithm.", 

599 "unit": "", 

600 }, 

601 ), 

602 ) 

603 pa_schema = pa_schema.append( 

604 pa.field( 

605 f"{b}_{alg_name}Flux", 

606 pa.float32(), 

607 nullable=b not in config.metadetect.shear_bands, 

608 metadata={ 

609 "doc": f"Flux in {b} band (measured with {alg_name} algorithm).", 

610 "unit": "", 

611 }, 

612 ), 

613 ) 

614 pa_schema = pa_schema.append( 

615 pa.field( 

616 f"{b}_{alg_name}FluxErr", 

617 pa.float32(), 

618 nullable=b not in config.metadetect.shear_bands, 

619 metadata={ 

620 "doc": f"Flux uncertainty in {b} band (measured with {alg_name} algorithm).", 

621 "unit": "", 

622 }, 

623 ), 

624 ) 

625 

626 return pa_schema 

627 

628 def runQuantum( 

629 self, 

630 qc: QuantumContext, 

631 inputRefs: InputQuantizedConnection, 

632 outputRefs: OutputQuantizedConnection, 

633 ) -> None: 

634 # Docstring inherited. 

635 

636 id_generator = self.config.id_generator.apply(qc.quantum.dataId) 

637 

638 if self.config.do_mask_bright_objects: 

639 ref_loader = ReferenceObjectLoader( 

640 dataIds=[ref.datasetRef.dataId for ref in inputRefs.ref_cat], 

641 refCats=[qc.get(ref) for ref in inputRefs.ref_cat], 

642 name=self.config.connections.ref_cat, 

643 config=self.config.ref_loader, 

644 log=self.log, 

645 ) 

646 ref_cat = ref_loader.loadRegion( 

647 qc.quantum.dataId.region, filterName=self.config.ref_loader_filter_name 

648 ) 

649 else: 

650 ref_cat = None 

651 

652 # Read the coadds and put them in the order defined by 

653 # config.photometry_bands (note that each MultipleCellCoadd object also 

654 # knows its own band, if that's needed). 

655 

656 coadds_by_band = { 

657 ref.dataId["band"]: qc.get(ref) 

658 for ref in inputRefs.input_coadds 

659 if ref.dataId["band"] in self.config.photometry_bands 

660 } 

661 

662 outputs = self.run( 

663 patch_coadds=coadds_by_band, 

664 id_generator=id_generator, 

665 ref_cat=ref_cat, 

666 ) 

667 qc.put(outputs, outputRefs) 

668 

669 def run( 

670 self, 

671 *, 

672 patch_coadds: Mapping[str, MultipleCellCoadd], 

673 id_generator: FullIdGenerator, 

674 ref_cat: SimpleCatalog | None, 

675 ) -> Struct: 

676 """Run metadetection on a patch. 

677 

678 Parameters 

679 ---------- 

680 patch_coadds : `~collections.abc.Mapping` [ \ 

681 `~lsst.cell_coadds.MultipleCellCoadd` ] 

682 Per-band, per-patch coadds, in the order specified by 

683 `MetadetectionShearConfig.photometry_bands`. 

684 id_generator : `~lsst.meas.base.FullIdGenerator` 

685 Generator for object IDs and to seed the random number generator. 

686 ref_cat : `lsst.afw.table.SimpleCatalog`, optional 

687 Reference catalog to use when masking bright stars. 

688 

689 Returns 

690 ------- 

691 results : `lsst.pipe.base.Struct` 

692 Structure with the following attributes: 

693 

694 - ``metadetect_catalog`` [ `pyarrow.Table` ]: the output object 

695 catalog for the patch, with schema equal to `metadetect_schema`. 

696 """ 

697 seed = id_generator.catalog_id 

698 self.rng = np.random.RandomState(seed) 

699 idstart = 0 

700 

701 grid = patch_coadds[self.config.photometry_bands[0]].grid 

702 nx_cells, ny_cells = grid.shape 

703 single_cell_tables: list[pa.Table] = [] 

704 for nx, ny in product(range(nx_cells), range(ny_cells)): 

705 cell_id = Index2D(nx, ny) 

706 bbox = grid.bbox_of(cell_id) 

707 cell_coadds = [patch_coadd.stitch(bbox) for patch_coadd in patch_coadds.values()] 

708 self.log.debug("Processing cell %s %s", nx, ny) 

709 

710 try: 

711 res = self.process_cell(cell_coadds, cell_id=cell_id) 

712 except Exception as e: 

713 self.log.error("Failed to process cell %s %s: %s", nx, ny, e) 

714 continue 

715 

716 if len(res) > 0: 

717 res["id"] = id_generator.arange(idstart, idstart + len(res)) 

718 # TODO: Avoid back and forth conversion between array and dict. 

719 da = self._dictify( 

720 res, 

721 tract=id_generator.data_id.tract.id, 

722 patch=id_generator.data_id.patch.id, 

723 ) 

724 table = pa.Table.from_pydict(da, self.metadetect_schema) 

725 

726 single_cell_tables.append(table) 

727 idstart += len(res) 

728 

729 if not single_cell_tables: 

730 raise MetadetectionProcessingError("No objects found in any cell") 

731 

732 # TODO: DM-53796 De-duplicate objects before concatenation. 

733 return Struct( 

734 metadetect_catalog=pa.concat_tables(single_cell_tables), 

735 ) 

736 

737 def process_cell( 

738 self, 

739 cell_coadds: Sequence[StitchedCoadd], 

740 cell_id: Index2D, 

741 ) -> pa.Table: 

742 """Run metadetection on a single cell. 

743 

744 Parameters 

745 ---------- 

746 cell_coadds : `~collections.abc.Sequence` [ \ 

747 `~lsst.cell_coadds.StitchedCoadd` ] 

748 Per-band, per-cell coadds, in the order specified by 

749 `MetadetectionShearConfig.photometry_bands`. 

750 cell_id : `~lsst.skymap.Index2D` 

751 The cell ID for the cell being processed. 

752 

753 Returns 

754 ------- 

755 metadetect_catalog : `pyarrow.Table` 

756 Output object catalog for the cell, with schema equal to 

757 `metadetect_schema`. 

758 """ 

759 

760 coadd_data = self._cell_to_coadd_data(cell_coadds) 

761 # TODO get bright star etc. info as input 

762 bright_info = [] 

763 

764 apply_apodized_edge_masks_mbexp(**coadd_data) 

765 

766 if len(bright_info) > 0: 

767 apply_apodized_bright_masks_mbexp(bright_info=bright_info, **coadd_data) 

768 

769 mask_frac = _get_mask_frac( 

770 coadd_data["mfrac_mbexp"], 

771 trim_pixels=0, 

772 ) 

773 

774 res = self.metadetect.run(rng=self.rng, **coadd_data) 

775 

776 comb_res = _make_comb_data( 

777 cell_coadd=cell_coadds[0], 

778 res=res, 

779 mask_frac=mask_frac, 

780 bands=[cell_coadd.band for cell_coadd in cell_coadds], 

781 cell_id=cell_id, 

782 ) 

783 

784 return comb_res 

785 

786 @staticmethod 

787 def _cell_to_coadd_data(cell_coadds: Sequence[StitchedCoadd]): 

788 coadd_data_list = [] 

789 for cell_coadd in cell_coadds: 

790 coadd_data = {} 

791 coadd_data["coadd_exp"] = cell_coadd.asExposure() 

792 coadd_data["coadd_noise_exp"] = cell_coadd.asExposure(noise_index=0) 

793 coadd_data["coadd_mfrac_exp"] = ExposureF(coadd_data["coadd_exp"], deep=True) 

794 coadd_data["coadd_mfrac_exp"].image = cell_coadd.mask_fractions 

795 coadd_data_list.append(coadd_data) 

796 

797 return extract_multiband_coadd_data(coadd_data_list) 

798 

799 def _dictify(self, data, tract: int, patch: int): 

800 output = {} 

801 # TODO: Move this to a better location after DP2. 

802 mapping = { 

803 "bmask": "bmask", 

804 "cell_x": "cell_x", 

805 "cell_y": "cell_y", 

806 "col": "x", 

807 "col_diff": "x_offset", # dropped. 

808 "dec": "coord_dec", 

809 "gauss_flags": "gauss_flags", 

810 "gauss_g_1": "gauss_g1", 

811 "gauss_g_2": "gauss_g2", 

812 "gauss_g_cov_11": "gauss_g1_g1_Cov", 

813 "gauss_g_cov_12": "gauss_g1_g2_Cov", # same as 21. 

814 "gauss_g_cov_22": "gauss_g2_g2_Cov", 

815 "gauss_obj_flags": "gauss_object_flags", 

816 "gauss_psf_flags": "gauss_psfReconvolved_flags", 

817 "gauss_psf_g_1": "gauss_psfReconvolved_g1", 

818 "gauss_psf_g_2": "gauss_psfReconvolved_g2", 

819 "gauss_psf_T": "gauss_psfReconvolved_T", 

820 "gauss_s2n": "gauss_snr", 

821 "gauss_T": "gauss_T", 

822 "gauss_T_err": "gauss_TErr", 

823 "gauss_T_flags": "gauss_T_flags", 

824 "gauss_T_ratio": "gauss_T_ratio", # dropped. 

825 "id": "shearObjectId", 

826 "mfrac": "mfrac", 

827 "ormask": "ormask", 

828 "pgauss_flags": "pgauss_flags", 

829 "pgauss_obj_flags": "pgauss_object_flags", 

830 "pgauss_s2n": "pgauss_snr", 

831 "pgauss_T": "pgauss_T", 

832 "pgauss_T_err": "pgauss_TErr", 

833 "pgauss_T_flags": "pgauss_T_flags", 

834 "pgauss_T_ratio": "pgauss_T_ratio", # dropped. 

835 "psfrec_flags": "psfOriginal_flags", 

836 "psfrec_g_1": "psfOriginal_e1", 

837 "psfrec_g_2": "psfOriginal_e2", 

838 "psfrec_T": "psfOriginal_T", 

839 "ra": "coord_ra", 

840 "row": "y", 

841 "row_diff": "y_offset", # dropped. 

842 "shear_type": "metaStep", 

843 "stamp_flags": "stamp_flags", 

844 } 

845 

846 for b, alg_name in product(self.config.photometry_bands, ("gauss", "pgauss")): 

847 mapping[f"{alg_name}_band_flux_{b}"] = f"{b}_{alg_name}Flux" 

848 mapping[f"{alg_name}_band_flux_err_{b}"] = f"{b}_{alg_name}FluxErr" 

849 mapping[f"{alg_name}_band_flux_flags_{b}"] = f"{b}_{alg_name}Flux_flags" 

850 

851 for name in mapping: 

852 if name in data.dtype.names: 

853 output[mapping.get(name, name)] = data[name] 

854 else: 

855 if "flags" in name.lower(): 

856 output[mapping.get(name, name)] = np.ones_like(data["id"], dtype=np.int32) 

857 else: 

858 output[mapping.get(name, name)] = np.ones_like(data["id"], dtype=np.float32) 

859 output[mapping.get(name, name)] *= np.nan 

860 

861 output["tract"] = tract * np.ones_like(data["id"], dtype=np.uint64) 

862 output["patch"] = patch * np.ones_like(data["id"], dtype=np.uint32) 

863 

864 return output 

865 

866 

867def _make_comb_data( 

868 cell_coadd, 

869 res, 

870 mask_frac, 

871 bands, 

872 cell_id, 

873): 

874 idinfo = cell_coadd.identifiers 

875 

876 copy_dt = [ 

877 # we will copy out of arrays to these 

878 ("psfrec_g_1", "f4"), 

879 ("psfrec_g_2", "f4"), 

880 ("gauss_psf_g_1", "f4"), 

881 ("gauss_psf_g_2", "f4"), 

882 ("gauss_g_1", "f4"), 

883 ("gauss_g_2", "f4"),