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

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.geom as geom 

43import lsst.pipe.base.connectionTypes as cT 

44from lsst.afw.image import ExposureF 

45from lsst.afw.table import SimpleCatalog 

46from lsst.cell_coadds import MultipleCellCoadd, StitchedCoadd 

47from lsst.daf.butler import DataCoordinate, DatasetRef 

48from lsst.meas.algorithms import LoadReferenceObjectsConfig, ReferenceObjectLoader 

49from lsst.meas.base import FullIdGenerator, SkyMapIdGeneratorConfig 

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

51from lsst.pipe.base import ( 

52 AlgorithmError, 

53 AnnotatedPartialOutputsError, 

54 InputQuantizedConnection, 

55 InvalidQuantumError, 

56 NoWorkFound, 

57 OutputQuantizedConnection, 

58 PipelineTask, 

59 PipelineTaskConfig, 

60 PipelineTaskConnections, 

61 QuantumContext, 

62 Struct, 

63) 

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

65from lsst.pipe.tasks.coaddBase import makeSkyInfo 

66from lsst.skymap import BaseSkyMap, Index2D 

67 

68 

69class MetadetectionProcessingError(AlgorithmError): 

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

71 

72 @property 

73 def metadata(self) -> dict: 

74 return {} 

75 

76 

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

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

79 

80 input_coadds = cT.Input( 

81 "deep_coadd_cell_predetection", 

82 storageClass="MultipleCellCoadd", 

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

84 multiple=True, 

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

86 ) 

87 

88 sky_map = cT.Input( 

89 doc="Cell-based skymap defining the patch structure.", 

90 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

91 storageClass="SkyMap", 

92 dimensions=("skymap",), 

93 ) 

94 

95 ref_cat = cT.PrerequisiteInput( 

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

97 name="the_monster_20250219", 

98 storageClass="SimpleCatalog", 

99 dimensions=("skypix",), 

100 deferLoad=True, 

101 multiple=True, 

102 ) 

103 

104 metadetect_catalog = cT.Output( 

105 "object_shear_patch", 

106 storageClass="ArrowTable", 

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

108 multiple=False, 

109 dimensions={"patch"}, 

110 ) 

111 

112 metadetect_schema = cT.InitOutput( 

113 "object_shear_schema", 

114 storageClass="ArrowSchema", 

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

116 ) 

117 

118 config: MetadetectionShearConfig 

119 

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

121 super().__init__(config=config) 

122 

123 if not config: 

124 return None 

125 

126 if not config.do_mask_bright_objects: 

127 del self.ref_cat 

128 

129 def adjustQuantum( 

130 self, 

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

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

133 label: str, 

134 data_id: DataCoordinate, 

135 ) -> tuple[ 

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

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

138 ]: 

139 # Docstring inherited. 

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

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

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

143 connection, original_input_coadds = inputs["input_coadds"] 

144 bands_missing = set(self.config.photometry_bands) 

145 adjusted_input_coadds = [] 

146 for ref in original_input_coadds: 

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

148 adjusted_input_coadds.append(ref) 

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

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

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

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

153 inputs.update(adjusted_inputs) 

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

155 return adjusted_inputs, {} 

156 

157 

158class MetadetectionShearConfig(PipelineTaskConfig, pipelineConnections=MetadetectionShearConnections): 

159 """Configuration definition for MetadetectionShearTask.""" 

160 

161 metadetect = ConfigurableField( 

162 target=MetadetectTask, 

163 doc="Configuration for metadetection.", 

164 ) 

165 

166 photometry_bands = ListField[str]( 

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

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

169 "not be processed.", 

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

171 ) 

172 

173 do_mask_bright_objects = Field[bool]( 

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

175 default=False, 

176 ) 

177 

178 ref_loader = ConfigField( 

179 dtype=LoadReferenceObjectsConfig, 

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

181 ) 

182 

183 ref_loader_filter_name = Field[str]( 

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

185 default="monster_DES_r", 

186 ) 

187 

188 border = Field[int]( 

189 "Border to apply to single cell images, if skymap has no cell borders", 

190 default=50, 

191 ) 

192 

193 do_cull_peaks_in_cell_borders = Field[bool]( 

194 "Cull detection peaks in the overlapping cell border regions?", 

195 default=True, 

196 ) 

197 

198 id_generator = SkyMapIdGeneratorConfig.make_field() 

199 

200 def setDefaults(self): 

201 super().setDefaults() 

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

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

204 

205 def validate(self): 

206 super().validate() 

207 if (shear_bands := self.metadetect.shear_bands) is not None and not set(shear_bands).issubset( 207 ↛ 210line 207 didn't jump to line 210 because the condition on line 207 was never true

208 self.photometry_bands 

209 ): 

210 raise FieldValidationError( 

211 self.__class__.metadetect, 

212 self, 

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

214 ) 

215 

216 

217class MetadetectionShearTask(PipelineTask): 

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

219 

220 _DefaultName: ClassVar[str] = "metadetectionShear" 

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

222 

223 config: MetadetectionShearConfig 

224 

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

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

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

228 self.makeSubtask("metadetect") 

229 

230 @classmethod 

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

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

233 

234 Parameters 

235 ---------- 

236 config : `MetadetectionShearConfig` 

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

238 

239 Returns 

240 ------- 

241 object_schema : `pyarrow.Schema` 

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

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

244 """ 

245 pa_schema = pa.schema( 

246 [ 

247 # Fields from pipeline bookkeeping. 

248 pa.field( 

249 "shearObjectId", 

250 pa.int64(), 

251 nullable=False, 

252 metadata={ 

253 "doc": ( 

254 "Unique identifier for a ShearObject, specific " 

255 "to a single metacalibration counterfactual image." 

256 ), 

257 "unit": "", 

258 }, 

259 ), 

260 pa.field( 

261 "tract", 

262 pa.int64(), 

263 nullable=False, 

264 metadata={ 

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

266 "unit": "", 

267 }, 

268 ), 

269 pa.field( 

270 "patch", 

271 pa.int64(), 

272 nullable=False, 

273 metadata={ 

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

275 "unit": "", 

276 }, 

277 ), 

278 pa.field( 

279 "cell_x", 

280 pa.int32(), 

281 nullable=False, 

282 metadata={ 

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

284 "unit": "", 

285 }, 

286 ), 

287 pa.field( 

288 "cell_y", 

289 pa.int32(), 

290 nullable=False, 

291 metadata={ 

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

293 "unit": "", 

294 }, 

295 ), 

296 pa.field( 

297 "is_cell_inner", 

298 pa.bool_(), 

299 nullable=False, 

300 metadata={ 

301 "doc": "Whether this object is within the inner region of the cell.", 

302 "unit": "", 

303 }, 

304 ), 

305 pa.field( 

306 "is_patch_inner", 

307 pa.bool_(), 

308 nullable=False, 

309 metadata={ 

310 "doc": "Whether this object is within the inner region of the patch.", 

311 "unit": "", 

312 }, 

313 ), 

314 pa.field( 

315 "is_tract_inner", 

316 pa.bool_(), 

317 nullable=False, 

318 metadata={ 

319 "doc": "Whether this object is within the inner region of the tract.", 

320 "unit": "", 

321 }, 

322 ), 

323 pa.field( 

324 "is_primary", 

325 pa.bool_(), 

326 nullable=False, 

327 metadata={ 

328 "doc": "Whether this object is to be selected after de-duplication", 

329 "unit": "", 

330 }, 

331 ), 

332 # Fields from metadetection (generic). 

333 pa.field( 

334 "metaStep", 

335 pa.string(), 

336 nullable=False, 

337 metadata={ 

338 "doc": ( 

339 "Type of artificial shear applied to image. " 

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

341 ), 

342 "unit": "", 

343 }, 

344 ), 

345 pa.field( 

346 "image_flags", 

347 pa.int32(), 

348 nullable=False, 

349 metadata={ 

350 "doc": "Flags for the image on which this measurement was made.", 

351 "unit": "", 

352 }, 

353 ), 

354 pa.field( 

355 "x", 

356 pa.float32(), 

357 nullable=False, 

358 metadata={ 

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

360 "unit": "", 

361 }, 

362 ), 

363 pa.field( 

364 "y", 

365 pa.float32(), 

366 nullable=False, 

367 metadata={ 

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

369 "unit": "", 

370 }, 

371 ), 

372 pa.field( 

373 "ra", 

374 pa.float64(), 

375 nullable=False, 

376 metadata={ 

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

378 "unit": "degrees", 

379 }, 

380 ), 

381 pa.field( 

382 "dec", 

383 pa.float64(), 

384 nullable=False, 

385 metadata={ 

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

387 "unit": "degrees", 

388 }, 

389 ), 

390 # Original PSF measurements 

391 pa.field( 

392 "psfOriginal_flags", 

393 pa.int32(), 

394 nullable=False, 

395 metadata={ 

396 "doc": "Flags for the original PSF measurement, ORed over the shear bands", 

397 "unit": "", 

398 }, 

399 ), 

400 pa.field( 

401 "psfOriginal_g1", 

402 pa.float32(), 

403 nullable=False, 

404 metadata={ 

405 "doc": ( 

406 "Reduced-shear g1 of the original PSF from adaptive moments, " 

407 "averaged over the shear bands." 

408 ), 

409 "unit": "", 

410 }, 

411 ), 

412 pa.field( 

413 "psfOriginal_g2", 

414 pa.float32(), 

415 nullable=False, 

416 metadata={ 

417 "doc": ( 

418 "Reduced-shear g2 of the original PSF from adaptive moments, " 

419 "averaged over the shear bands." 

420 ), 

421 "unit": "", 

422 }, 

423 ), 

424 pa.field( 

425 "psfOriginal_T", 

426 pa.float32(), 

427 nullable=False, 

428 metadata={ 

429 "doc": ( 

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

431 "averaged over the shear bands." 

432 ), 

433 "unit": "arcseconds squared", 

434 }, 

435 ), 

436 pa.field( 

437 "bmask_flags", 

438 pa.int32(), 

439 nullable=False, 

440 metadata={ 

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

442 "unit": "", 

443 }, 

444 ), 

445 pa.field( 

446 "ormask_flags", 

447 pa.int32(), 

448 nullable=False, 

449 metadata={ 

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

451 "unit": "", 

452 }, 

453 ), 

454 pa.field( 

455 "mfrac", 

456 pa.float32(), 

457 nullable=False, 

458 metadata={ 

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

460 "unit": "", 

461 }, 

462 ), 

463 # Fields that come only from gauss algorithm. 

464 # Reconvolved PSF measurements (gauss) 

465 pa.field( 

466 "gauss_psfReconvolved_flags", 

467 pa.int32(), 

468 nullable=False, 

469 metadata={ 

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

471 "unit": "", 

472 }, 

473 ), 

474 pa.field( 

475 "gauss_psfReconvolved_g1", 

476 pa.float32(), 

477 nullable=False, 

478 metadata={ 

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

480 "unit": "", 

481 }, 

482 ), 

483 pa.field( 

484 "gauss_psfReconvolved_g2", 

485 pa.float32(), 

486 nullable=False, 

487 metadata={ 

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

489 "unit": "", 

490 }, 

491 ), 

492 pa.field( 

493 "gauss_psfReconvolved_T", 

494 pa.float32(), 

495 nullable=False, 

496 metadata={ 

497 "doc": ( 

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

499 ), 

500 "unit": "arcseconds squared", 

501 }, 

502 ), 

503 # Object measurements (gauss algorithm). 

504 pa.field( 

505 "gauss_g1", 

506 pa.float32(), 

507 nullable=False, 

508 metadata={ 

509 "doc": ( 

510 "Reduced-shear g1 measurement of the ShearObject " 

511 "(measured with gauss algorithm)." 

512 ), 

513 "unit": "", 

514 }, 

515 ), 

516 pa.field( 

517 "gauss_g2", 

518 pa.float32(), 

519 nullable=False, 

520 metadata={ 

521 "doc": ( 

522 "Reduced-shear g2 measurement of the ShearObject " 

523 "(measured with gauss algorithm)." 

524 ), 

525 "unit": "", 

526 }, 

527 ), 

528 pa.field( 

529 "gauss_g1_g1_Cov", 

530 pa.float32(), 

531 nullable=False, 

532 metadata={ 

533 "doc": ( 

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

535 "(measured with gauss algorithm)." 

536 ), 

537 "unit": "", 

538 }, 

539 ), 

540 pa.field( 

541 "gauss_g1_g2_Cov", 

542 pa.float32(), 

543 nullable=False, 

544 metadata={ 

545 "doc": ( 

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

547 "(measured with gauss algorithm)." 

548 ), 

549 "unit": "", 

550 }, 

551 ), 

552 pa.field( 

553 "gauss_g2_g2_Cov", 

554 pa.float32(), 

555 nullable=False, 

556 metadata={ 

557 "doc": ( 

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

559 "(measured with gauss algorithm)." 

560 ), 

561 "unit": "", 

562 }, 

563 ), 

564 ], 

565 metadata={ 

566 "shear_step": str(SHEAR_STEP), 

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

568 }, 

569 ) 

570 

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

572 pa_schema = pa_schema.append( 

573 pa.field( 

574 f"{alg_name}_snr", 

575 pa.float32(), 

576 nullable=False, 

577 metadata={ 

578 "doc": ( 

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

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

581 ), 

582 "unit": "", 

583 }, 

584 ), 

585 ) 

586 pa_schema = pa_schema.append( 

587 pa.field( 

588 f"{alg_name}_T", 

589 pa.float32(), 

590 nullable=False, 

591 metadata={ 

592 "doc": ( 

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

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

595 ), 

596 "unit": "arcseconds squared", 

597 }, 

598 ), 

599 ) 

600 pa_schema = pa_schema.append( 

601 pa.field( 

602 f"{alg_name}_TErr", 

603 pa.float32(), 

604 nullable=False, 

605 metadata={ 

606 "doc": ( 

607 "Uncertainty in the trace measurement of the ShearObject " 

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

609 ), 

610 "unit": "arcseconds squared", 

611 }, 

612 ), 

613 ) 

614 pa_schema = pa_schema.append( 

615 pa.field( 

616 f"{alg_name}_shape_flags", 

617 pa.int32(), 

618 nullable=False, 

619 metadata={ 

620 "doc": ( 

621 "Flags for the second order moments measurement of the ShearObject " 

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

623 ), 

624 "unit": "", 

625 }, 

626 ), 

627 ) 

628 pa_schema = pa_schema.append( 

629 pa.field( 

630 f"{alg_name}_object_flags", 

631 pa.int32(), 

632 nullable=False, 

633 metadata={ 

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

635 "unit": "", 

636 }, 

637 ), 

638 ) 

639 pa_schema = pa_schema.append( 

640 pa.field( 

641 f"{alg_name}_flags", 

642 pa.int32(), 

643 nullable=False, 

644 metadata={ 

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

646 "unit": "", 

647 }, 

648 ), 

649 ) 

650 

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

652 for b in config.photometry_bands: 

653 pa_schema = pa_schema.append( 

654 pa.field( 

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

656 pa.int32(), 

657 nullable=False, 

658 metadata={ 

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

660 "unit": "", 

661 }, 

662 ), 

663 ) 

664 pa_schema = pa_schema.append( 

665 pa.field( 

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

667 pa.float32(), 

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

669 metadata={ 

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

671 "unit": "", 

672 }, 

673 ), 

674 ) 

675 pa_schema = pa_schema.append( 

676 pa.field( 

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

678 pa.float32(), 

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

680 metadata={ 

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

682 "unit": "", 

683 }, 

684 ), 

685 ) 

686 

687 return pa_schema 

688 

689 def validate_skymap_config(self, skymap_config: BaseSkyMap.ConfigClass) -> None: 

690 if not skymap_config.tractBuilder.name == "cells": 

691 raise InvalidQuantumError("MetadetectionShearTask requires a cell-based skymap.") 

692 

693 cell_config = skymap_config.tractBuilder.active 

694 if (self.config.border == 0 and cell_config.cellBorder == 0) or ( 

695 self.config.border > 0 and cell_config.cellBorder > 0 

696 ): 

697 raise InvalidQuantumError( 

698 "MetadetectionShearTask requires a positive border to be set either in the skymap config " 

699 "or in the task config (but not in both)." 

700 ) 

701 

702 if self.config.border: 

703 # In case cellInnerDimensions are different in different directions 

704 # (which is rare in practice), take the min of it 

705 cell_inner_dimensions = min(cell_config.cellInnerDimensions) 

706 if self.config.border > ( 

707 max_border_value := cell_config.numCellsInPatchBorder * cell_inner_dimensions 

708 ): 

709 raise InvalidQuantumError( 

710 "The border value is too large for the skymap configuration. " 

711 f"Maximum border value is {max_border_value}." 

712 ) 

713 

714 # Ensure that the amount by which we withdraw inwards does not 

715 # create gaps at tract boundaries. 

716 # tractOverlap is specified in degrees. 

717 if ( 

718 skymap_config.tractOverlap * 3600 / skymap_config.pixelScale 

719 < self.count_cells_along_edges(skymap_config) * cell_inner_dimensions + cell_config.cellBorder 

720 ): 

721 raise InvalidQuantumError( 

722 "The tract overlap is insufficient given the borders. " 

723 "This will result in missed regions between adjacent tracts." 

724 ) 

725 

726 @staticmethod 

727 def count_cells_along_edges(skymap_config: BaseSkyMap.ConfigClass) -> int: 

728 """Count the number of cells along the edges to skip processing.""" 

729 if skymap_config.tractBuilder["cells"].cellBorder > 0: 

730 return 0 

731 

732 return skymap_config.tractBuilder["cells"].numCellsInPatchBorder 

733 

734 def runQuantum( 

735 self, 

736 qc: QuantumContext, 

737 inputRefs: InputQuantizedConnection, 

738 outputRefs: OutputQuantizedConnection, 

739 ) -> None: 

740 # Docstring inherited. 

741 

742 # Get the skyMap for this quantum 

743 sky_map = qc.get(inputRefs.sky_map) 

744 

745 self.validate_skymap_config(sky_map.config) 

746 

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

748 

749 if self.config.do_mask_bright_objects: 

750 ref_loader = ReferenceObjectLoader( 

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

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

753 name=self.config.connections.ref_cat, 

754 config=self.config.ref_loader, 

755 log=self.log, 

756 ) 

757 ref_cat = ref_loader.loadRegion( 

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

759 ) 

760 else: 

761 ref_cat = None 

762 

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

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

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

766 

767 coadds_by_band = { 

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

769 for ref in inputRefs.input_coadds 

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

771 } 

772 

773 try: 

774 outputs = self.run( 

775 patch_coadds=coadds_by_band, 

776 id_generator=id_generator, 

777 sky_map=sky_map, 

778 ref_cat=ref_cat, 

779 ) 

780 except AlgorithmError as err: 

781 # We know there are no actual outputs in this case, but this is 

782 # still the right exception to raise (it's just badly named). 

783 raise AnnotatedPartialOutputsError.annotate(err, self, log=self.log) from err 

784 qc.put(outputs, outputRefs) 

785 

786 def run( 

787 self, 

788 *, 

789 patch_coadds: Mapping[str, MultipleCellCoadd], 

790 id_generator: FullIdGenerator, 

791 sky_map: BaseSkyMap, 

792 ref_cat: SimpleCatalog | None, 

793 ) -> Struct: 

794 """Run metadetection on a patch. 

795 

796 Parameters 

797 ---------- 

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

799 `~lsst.cell_coadds.MultipleCellCoadd` ] 

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

801 `MetadetectionShearConfig.photometry_bands`. 

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

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

804 sky_map : `~lsst.skymap.BaseSkyMap` 

805 Sky map to use for determining the patch boundaries. 

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

807 Reference catalog to use when masking bright stars. 

808 

809 Returns 

810 ------- 

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

812 Structure with the following attributes: 

813 

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

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

816 """ 

817 seed = id_generator.catalog_id 

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

819 idstart = 0 

820 

821 sky_info = makeSkyInfo( 

822 sky_map, 

823 tractId=id_generator.data_id.tract.id, 

824 patchId=id_generator.data_id.patch.id, 

825 ) 

826 # Since cells have no borders, we cannot apply the metacal 

827 # procedure to the edge cells in the same way as inner cells. 

828 # We can do so for all that cells that are marked as borders, 

829 # but it has to be at least 1. 

830 if (num_edge_cells_skip := sky_map.config.tractBuilder.active.numCellsInPatchBorder) < 1: 

831 raise InvalidQuantumError("No border cells found in the skymap configuration.") 

832 

833 dilate_by = self.config.border or sky_map.config.tractBuilder.active.cellBorder 

834 border = dilate_by if self.config.do_cull_peaks_in_cell_borders else 0 

835 

836 grid = patch_coadds[self.config.metadetect.shear_bands[0]].grid 

837 # Undo any padding that was applied when creating the grid. 

838 # We only support metadetection on non-overlapping single cells. 

839 # lsst_cells_v1 has overlapping cells, so we need to undo the padding. 

840 grid = grid.from_bbox_shape(bbox=grid.bbox, shape=grid.shape, padding=0) 

841 

842 nx_cells, ny_cells = grid.shape 

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

844 for nx, ny in product( 

845 range(num_edge_cells_skip, nx_cells - num_edge_cells_skip), 

846 range(num_edge_cells_skip, ny_cells - num_edge_cells_skip), 

847 ): 

848 cell_id = Index2D(nx, ny) 

849 bbox = grid.bbox_of(cell_id).dilatedBy(dilate_by) 

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

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

852 

853 try: 

854 res = self.process_cell(cell_coadds, cell_id=cell_id, border=border) 

855 except Exception as e: 

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

857 continue 

858 

859 if len(res) > 0: 

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

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

862 da = self._dictify( 

863 res, 

864 tract=id_generator.data_id.tract.id, 

865 patch=id_generator.data_id.patch.id, 

866 ) 

867 

868 self._fill_is_inner(da, sky_info.patchInfo[cell_id].inner_bbox, "is_cell_inner") 

869 self._fill_is_inner(da, sky_info.patchInfo.inner_bbox, "is_patch_inner") 

870 # Calculating is_tract_inner follows a different logic with 

871 # regions instead of boxes. 

872 da["is_tract_inner"] = sky_info.tractInfo.inner_sky_region.contains( 

873 da["ra"] * geom.degrees, 

874 da["dec"] * geom.degrees, 

875 ) 

876 da["is_primary"] = da["is_cell_inner"] & da["is_patch_inner"] & da["is_tract_inner"] 

877 

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

879 

880 single_cell_tables.append(table) 

881 idstart += len(res) 

882 

883 if not single_cell_tables: 

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

885 

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

887 return Struct( 

888 metadetect_catalog=pa.concat_tables(single_cell_tables), 

889 ) 

890 

891 def process_cell( 

892 self, 

893 cell_coadds: Sequence[StitchedCoadd],