Coverage for python/lsst/pipe/tasks/finalizeCharacterization.py: 15%

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"""Task to run a finalized image characterization, using additional data. 

23""" 

24 

25__all__ = [ 

26 'FinalizeCharacterizationConnections', 

27 'FinalizeCharacterizationConfig', 

28 'FinalizeCharacterizationTask', 

29 'FinalizeCharacterizationDetectorConnections', 

30 'FinalizeCharacterizationDetectorConfig', 

31 'FinalizeCharacterizationDetectorTask', 

32 'ConsolidateFinalizeCharacterizationDetectorConnections', 

33 'ConsolidateFinalizeCharacterizationDetectorConfig', 

34 'ConsolidateFinalizeCharacterizationDetectorTask', 

35] 

36 

37import astropy.table 

38import astropy.units as u 

39import numpy as np 

40import esutil 

41from smatch.matcher import Matcher 

42 

43 

44import lsst.pex.config as pexConfig 

45import lsst.pipe.base as pipeBase 

46import lsst.daf.base as dafBase 

47import lsst.afw.table as afwTable 

48import lsst.meas.algorithms as measAlg 

49import lsst.meas.extensions.piff.piffPsfDeterminer # noqa: F401 

50from lsst.meas.algorithms import MeasureApCorrTask 

51from lsst.meas.base import SingleFrameMeasurementTask, ApplyApCorrTask 

52from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry 

53 

54from .reserveIsolatedStars import ReserveIsolatedStarsTask 

55from lsst.obs.base.utils import TableVStack 

56 

57 

58class FinalizeCharacterizationConnectionsBase( 

59 pipeBase.PipelineTaskConnections, 

60 dimensions=('instrument', 'visit',), 

61 defaultTemplates={}, 

62): 

63 src_schema = pipeBase.connectionTypes.InitInput( 

64 doc='Input schema used for src catalogs.', 

65 name='src_schema', 

66 storageClass='SourceCatalog', 

67 ) 

68 isolated_star_cats = pipeBase.connectionTypes.Input( 

69 doc=('Catalog of isolated stars with average positions, number of associated ' 

70 'sources, and indexes to the isolated_star_sources catalogs.'), 

71 name='isolated_star_presource_associations', 

72 storageClass='ArrowAstropy', 

73 dimensions=('instrument', 'tract', 'skymap'), 

74 deferLoad=True, 

75 multiple=True, 

76 ) 

77 isolated_star_sources = pipeBase.connectionTypes.Input( 

78 doc=('Catalog of isolated star sources with sourceIds, and indexes to the ' 

79 'isolated_star_cats catalogs.'), 

80 name='isolated_star_presources', 

81 storageClass='ArrowAstropy', 

82 dimensions=('instrument', 'tract', 'skymap'), 

83 deferLoad=True, 

84 multiple=True, 

85 ) 

86 fgcm_standard_star = pipeBase.connectionTypes.Input( 

87 doc=('Catalog of fgcm for color corrections, and indexes to the ' 

88 'isolated_star_cats catalogs.'), 

89 name='fgcm_standard_star', 

90 storageClass='ArrowAstropy', 

91 dimensions=('instrument', 'tract', 'skymap'), 

92 deferLoad=True, 

93 multiple=True, 

94 ) 

95 

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

97 super().__init__(config=config) 

98 if config is None: 

99 return None 

100 if not self.config.psf_determiner['piff'].useColor: 

101 self.inputs.remove("fgcm_standard_star") 

102 

103 

104class FinalizeCharacterizationConnections( 

105 FinalizeCharacterizationConnectionsBase, 

106 dimensions=('instrument', 'visit',), 

107 defaultTemplates={}, 

108): 

109 srcs = pipeBase.connectionTypes.Input( 

110 doc='Source catalogs for the visit', 

111 name='src', 

112 storageClass='SourceCatalog', 

113 dimensions=('instrument', 'visit', 'detector'), 

114 deferLoad=True, 

115 multiple=True, 

116 deferGraphConstraint=True, 

117 ) 

118 calexps = pipeBase.connectionTypes.Input( 

119 doc='Calexps for the visit', 

120 name='calexp', 

121 storageClass='ExposureF', 

122 dimensions=('instrument', 'visit', 'detector'), 

123 deferLoad=True, 

124 multiple=True, 

125 ) 

126 finalized_psf_ap_corr_cat = pipeBase.connectionTypes.Output( 

127 doc=('Per-visit finalized psf models and aperture corrections. This ' 

128 'catalog uses detector id for the id and are sorted for fast ' 

129 'lookups of a detector.'), 

130 name='finalized_psf_ap_corr_catalog', 

131 storageClass='ExposureCatalog', 

132 dimensions=('instrument', 'visit'), 

133 ) 

134 finalized_src_table = pipeBase.connectionTypes.Output( 

135 doc=('Per-visit catalog of measurements for psf/flag/etc.'), 

136 name='finalized_src_table', 

137 storageClass='ArrowAstropy', 

138 dimensions=('instrument', 'visit'), 

139 ) 

140 

141 

142class FinalizeCharacterizationDetectorConnections( 

143 FinalizeCharacterizationConnectionsBase, 

144 dimensions=('instrument', 'visit', 'detector',), 

145 defaultTemplates={}, 

146): 

147 src = pipeBase.connectionTypes.Input( 

148 doc='Source catalog for the visit/detector.', 

149 name='src', 

150 storageClass='SourceCatalog', 

151 dimensions=('instrument', 'visit', 'detector'), 

152 ) 

153 calexp = pipeBase.connectionTypes.Input( 

154 doc='Calibrated exposure for the visit/detector.', 

155 name='calexp', 

156 storageClass='ExposureF', 

157 dimensions=('instrument', 'visit', 'detector'), 

158 ) 

159 finalized_psf_ap_corr_detector_cat = pipeBase.connectionTypes.Output( 

160 doc=('Per-visit/per-detector finalized psf models and aperture corrections. This ' 

161 'catalog uses detector id for the id.'), 

162 name='finalized_psf_ap_corr_detector_catalog', 

163 storageClass='ExposureCatalog', 

164 dimensions=('instrument', 'visit', 'detector'), 

165 ) 

166 finalized_src_detector_table = pipeBase.connectionTypes.Output( 

167 doc=('Per-visit/per-detector catalog of measurements for psf/flag/etc.'), 

168 name='finalized_src_detector_table', 

169 storageClass='ArrowAstropy', 

170 dimensions=('instrument', 'visit', 'detector'), 

171 ) 

172 

173 

174class FinalizeCharacterizationConfigBase( 

175 pipeBase.PipelineTaskConfig, 

176 pipelineConnections=FinalizeCharacterizationConnectionsBase, 

177): 

178 """Configuration for FinalizeCharacterizationBaseTask.""" 

179 source_selector = sourceSelectorRegistry.makeField( 

180 doc="How to select sources", 

181 default="science" 

182 ) 

183 id_column = pexConfig.Field( 

184 doc='Name of column in isolated_star_sources with source id.', 

185 dtype=str, 

186 default='sourceId', 

187 ) 

188 reserve_selection = pexConfig.ConfigurableField( 

189 target=ReserveIsolatedStarsTask, 

190 doc='Task to select reserved stars', 

191 ) 

192 make_psf_candidates = pexConfig.ConfigurableField( 

193 target=measAlg.MakePsfCandidatesTask, 

194 doc='Task to make psf candidates from selected stars.', 

195 ) 

196 psf_determiner = measAlg.psfDeterminerRegistry.makeField( 

197 'PSF Determination algorithm', 

198 default='piff' 

199 ) 

200 measurement = pexConfig.ConfigurableField( 

201 target=SingleFrameMeasurementTask, 

202 doc='Measure sources for aperture corrections' 

203 ) 

204 measure_ap_corr = pexConfig.ConfigurableField( 

205 target=MeasureApCorrTask, 

206 doc="Subtask to measure aperture corrections" 

207 ) 

208 apply_ap_corr = pexConfig.ConfigurableField( 

209 target=ApplyApCorrTask, 

210 doc="Subtask to apply aperture corrections" 

211 ) 

212 

213 def setDefaults(self): 

214 super().setDefaults() 

215 

216 source_selector = self.source_selector['science'] 

217 source_selector.setDefaults() 

218 

219 # We use the source selector only to select out flagged objects 

220 # and signal-to-noise. Isolated, unresolved sources are handled 

221 # by the isolated star catalog. 

222 

223 source_selector.doFlags = True 

224 source_selector.doSignalToNoise = True 

225 source_selector.doFluxLimit = False 

226 source_selector.doUnresolved = False 

227 source_selector.doIsolated = False 

228 

229 source_selector.signalToNoise.minimum = 50.0 

230 source_selector.signalToNoise.maximum = 1000.0 

231 

232 source_selector.signalToNoise.fluxField = 'base_GaussianFlux_instFlux' 

233 source_selector.signalToNoise.errField = 'base_GaussianFlux_instFluxErr' 

234 

235 source_selector.flags.bad = ['base_PixelFlags_flag_edge', 

236 'base_PixelFlags_flag_nodata', 

237 'base_PixelFlags_flag_interpolatedCenter', 

238 'base_PixelFlags_flag_saturatedCenter', 

239 'base_PixelFlags_flag_crCenter', 

240 'base_PixelFlags_flag_bad', 

241 'base_PixelFlags_flag_interpolated', 

242 'base_PixelFlags_flag_saturated', 

243 'slot_Centroid_flag', 

244 'base_GaussianFlux_flag'] 

245 

246 # Configure aperture correction to select only high s/n sources (that 

247 # were used in the psf modeling) to avoid background problems when 

248 # computing the aperture correction map. 

249 self.measure_ap_corr.sourceSelector = 'science' 

250 

251 ap_selector = self.measure_ap_corr.sourceSelector['science'] 

252 # We do not need to filter flags or unresolved because we have used 

253 # the filtered isolated stars as an input 

254 ap_selector.doFlags = False 

255 ap_selector.doUnresolved = False 

256 

257 import lsst.meas.modelfit # noqa: F401 

258 import lsst.meas.extensions.photometryKron # noqa: F401 

259 import lsst.meas.extensions.convolved # noqa: F401 

260 import lsst.meas.extensions.gaap # noqa: F401 

261 import lsst.meas.extensions.shapeHSM # noqa: F401 

262 

263 # Set up measurement defaults 

264 self.measurement.plugins.names = [ 

265 'base_FPPosition', 

266 'base_PsfFlux', 

267 'base_GaussianFlux', 

268 'modelfit_DoubleShapeletPsfApprox', 

269 'modelfit_CModel', 

270 'ext_photometryKron_KronFlux', 

271 'ext_convolved_ConvolvedFlux', 

272 'ext_gaap_GaapFlux', 

273 'ext_shapeHSM_HsmShapeRegauss', 

274 'ext_shapeHSM_HsmSourceMoments', 

275 'ext_shapeHSM_HsmPsfMoments', 

276 'ext_shapeHSM_HsmSourceMomentsRound', 

277 'ext_shapeHSM_HigherOrderMomentsSource', 

278 'ext_shapeHSM_HigherOrderMomentsPSF', 

279 ] 

280 self.measurement.slots.modelFlux = 'modelfit_CModel' 

281 self.measurement.plugins['ext_convolved_ConvolvedFlux'].seeing.append(8.0) 

282 self.measurement.plugins['ext_gaap_GaapFlux'].sigmas = [ 

283 0.5, 

284 0.7, 

285 1.0, 

286 1.5, 

287 2.5, 

288 3.0 

289 ] 

290 self.measurement.plugins['ext_gaap_GaapFlux'].doPsfPhotometry = True 

291 self.measurement.slots.shape = 'ext_shapeHSM_HsmSourceMoments' 

292 self.measurement.slots.psfShape = 'ext_shapeHSM_HsmPsfMoments' 

293 self.measurement.plugins['ext_shapeHSM_HsmShapeRegauss'].deblendNChild = "" 

294 

295 # TODO: Remove in DM-44658, streak masking to happen only in ip_diffim 

296 # Keep track of which footprints contain streaks 

297 self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['STREAK'] 

298 self.measurement.plugins['base_PixelFlags'].masksFpCenter = ['STREAK'] 

299 

300 # Turn off slot setting for measurement for centroid and shape 

301 # (for which we use the input src catalog measurements) 

302 self.measurement.slots.centroid = None 

303 self.measurement.slots.apFlux = None 

304 self.measurement.slots.calibFlux = None 

305 

306 names = self.measurement.plugins['ext_convolved_ConvolvedFlux'].getAllResultNames() 

307 self.measure_ap_corr.allowFailure += names 

308 names = self.measurement.plugins["ext_gaap_GaapFlux"].getAllGaapResultNames() 

309 self.measure_ap_corr.allowFailure += names 

310 

311 

312class FinalizeCharacterizationConfig( 

313 FinalizeCharacterizationConfigBase, 

314 pipelineConnections=FinalizeCharacterizationConnections, 

315): 

316 pass 

317 

318 

319class FinalizeCharacterizationDetectorConfig( 

320 FinalizeCharacterizationConfigBase, 

321 pipelineConnections=FinalizeCharacterizationDetectorConnections, 

322): 

323 pass 

324 

325 

326class FinalizeCharacterizationTaskBase(pipeBase.PipelineTask): 

327 """Run final characterization on exposures.""" 

328 ConfigClass = FinalizeCharacterizationConfigBase 

329 _DefaultName = 'finalize_characterization_base' 

330 

331 def __init__(self, initInputs=None, **kwargs): 

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

333 

334 self.schema_mapper, self.schema = self._make_output_schema_mapper( 

335 initInputs['src_schema'].schema 

336 ) 

337 

338 self.makeSubtask('reserve_selection') 

339 self.makeSubtask('source_selector') 

340 self.makeSubtask('make_psf_candidates') 

341 self.makeSubtask('psf_determiner') 

342 self.makeSubtask('measurement', schema=self.schema) 

343 self.makeSubtask('measure_ap_corr', schema=self.schema) 

344 self.makeSubtask('apply_ap_corr', schema=self.schema) 

345 

346 # Only log warning and fatal errors from the source_selector 

347 self.source_selector.log.setLevel(self.source_selector.log.WARN) 

348 self.isPsfDeterminerPiff = False 

349 if isinstance(self.psf_determiner, lsst.meas.extensions.piff.piffPsfDeterminer.PiffPsfDeterminerTask): 

350 self.isPsfDeterminerPiff = True 

351 

352 def _make_output_schema_mapper(self, input_schema): 

353 """Make the schema mapper from the input schema to the output schema. 

354 

355 Parameters 

356 ---------- 

357 input_schema : `lsst.afw.table.Schema` 

358 Input schema. 

359 

360 Returns 

361 ------- 

362 mapper : `lsst.afw.table.SchemaMapper` 

363 Schema mapper 

364 output_schema : `lsst.afw.table.Schema` 

365 Output schema (with alias map) 

366 """ 

367 mapper = afwTable.SchemaMapper(input_schema) 

368 mapper.addMinimalSchema(afwTable.SourceTable.makeMinimalSchema()) 

369 mapper.addMapping(input_schema['slot_Centroid_x'].asKey()) 

370 mapper.addMapping(input_schema['slot_Centroid_y'].asKey()) 

371 

372 # The aperture fields may be used by the psf determiner. 

373 aper_fields = input_schema.extract('base_CircularApertureFlux_*') 

374 for field, item in aper_fields.items(): 

375 mapper.addMapping(item.key) 

376 

377 # The following two may be redundant, but then the mapping is a no-op. 

378 # Note that the slot_CalibFlux mapping will copy over any 

379 # normalized compensated fluxes that are used for calibration. 

380 apflux_fields = input_schema.extract('slot_ApFlux_*') 

381 for field, item in apflux_fields.items(): 

382 mapper.addMapping(item.key) 

383 

384 calibflux_fields = input_schema.extract('slot_CalibFlux_*') 

385 for field, item in calibflux_fields.items(): 

386 mapper.addMapping(item.key) 

387 

388 mapper.addMapping( 

389 input_schema[self.config.source_selector.active.signalToNoise.fluxField].asKey(), 

390 'calib_psf_selection_flux') 

391 mapper.addMapping( 

392 input_schema[self.config.source_selector.active.signalToNoise.errField].asKey(), 

393 'calib_psf_selection_flux_err') 

394 

395 output_schema = mapper.getOutputSchema() 

396 

397 output_schema.addField( 

398 'calib_psf_candidate', 

399 type='Flag', 

400 doc=('set if the source was a candidate for PSF determination, ' 

401 'as determined from FinalizeCharacterizationTask.'), 

402 ) 

403 output_schema.addField( 

404 'calib_psf_reserved', 

405 type='Flag', 

406 doc=('set if source was reserved from PSF determination by ' 

407 'FinalizeCharacterizationTask.'), 

408 ) 

409 output_schema.addField( 

410 'calib_psf_used', 

411 type='Flag', 

412 doc=('set if source was used in the PSF determination by ' 

413 'FinalizeCharacterizationTask.'), 

414 ) 

415 output_schema.addField( 

416 'visit', 

417 type=np.int64, 

418 doc='Visit number for the sources.', 

419 ) 

420 output_schema.addField( 

421 'detector', 

422 type=np.int32, 

423 doc='Detector number for the sources.', 

424 ) 

425 output_schema.addField( 

426 'psf_color_value', 

427 type=np.float32, 

428 doc="Color used in PSF fit." 

429 ) 

430 output_schema.addField( 

431 'psf_color_type', 

432 type=str, 

433 size=10, 

434 doc="Color used in PSF fit." 

435 ) 

436 output_schema.addField( 

437 'psf_max_value', 

438 type=np.float32, 

439 doc="Maximum value in the star image used to train PSF.", 

440 doReplace=True, 

441 ) 

442 

443 alias_map = input_schema.getAliasMap() 

444 alias_map_output = afwTable.AliasMap() 

445 alias_map_output.set('slot_Centroid', alias_map.get('slot_Centroid')) 

446 alias_map_output.set('slot_ApFlux', alias_map.get('slot_ApFlux')) 

447 alias_map_output.set('slot_CalibFlux', alias_map.get('slot_CalibFlux')) 

448 

449 output_schema.setAliasMap(alias_map_output) 

450 

451 return mapper, output_schema 

452 

453 def _make_selection_schema_mapper(self, input_schema): 

454 """Make the schema mapper from the input schema to the selection schema. 

455 

456 Parameters 

457 ---------- 

458 input_schema : `lsst.afw.table.Schema` 

459 Input schema. 

460 

461 Returns 

462 ------- 

463 mapper : `lsst.afw.table.SchemaMapper` 

464 Schema mapper 

465 selection_schema : `lsst.afw.table.Schema` 

466 Selection schema (with alias map) 

467 """ 

468 mapper = afwTable.SchemaMapper(input_schema) 

469 mapper.addMinimalSchema(input_schema) 

470 

471 selection_schema = mapper.getOutputSchema() 

472 

473 selection_schema.setAliasMap(input_schema.getAliasMap()) 

474 

475 selection_schema.addField( 

476 'psf_color_value', 

477 type=np.float32, 

478 doc="Color used in PSF fit." 

479 ) 

480 selection_schema.addField( 

481 'psf_color_type', 

482 type=str, 

483 size=10, 

484 doc="Color used in PSF fit." 

485 ) 

486 selection_schema.addField( 

487 'psf_max_value', 

488 type=np.float32, 

489 doc="Maximum value in the star image used to train PSF.", 

490 doReplace=True, 

491 ) 

492 

493 return mapper, selection_schema 

494 

495 def concat_isolated_star_cats(self, band, isolated_star_cat_dict, isolated_star_source_dict): 

496 """ 

497 Concatenate isolated star catalogs and make reserve selection. 

498 

499 Parameters 

500 ---------- 

501 band : `str` 

502 Band name. Used to select reserved stars. 

503 isolated_star_cat_dict : `dict` 

504 Per-tract dict of isolated star catalog handles. 

505 isolated_star_source_dict : `dict` 

506 Per-tract dict of isolated star source catalog handles. 

507 

508 Returns 

509 ------- 

510 isolated_table : `np.ndarray` (N,) 

511 Table of isolated stars, with indexes to isolated sources. 

512 Returns None if there are no usable isolated catalogs. 

513 isolated_source_table : `np.ndarray` (M,) 

514 Table of isolated sources, with indexes to isolated stars. 

515 Returns None if there are no usable isolated catalogs. 

516 """ 

517 isolated_tables = [] 

518 isolated_sources = [] 

519 merge_cat_counter = 0 

520 merge_source_counter = 0 

521 

522 for tract in isolated_star_cat_dict: 

523 astropy_cat = isolated_star_cat_dict[tract].get() 

524 table_cat = np.asarray(astropy_cat) 

525 

526 astropy_source = isolated_star_source_dict[tract].get( 

527 parameters={'columns': [self.config.id_column, 'obj_index']} 

528 ) 

529 table_source = np.asarray(astropy_source) 

530 

531 # Cut isolated star table to those observed in this band, and adjust indexes 

532 (use_band,) = (table_cat[f'nsource_{band}'] > 0).nonzero() 

533 

534 if len(use_band) == 0: 

535 # There are no sources in this band in this tract. 

536 self.log.info("No sources found in %s band in tract %d.", band, tract) 

537 continue 

538 

539 # With the following matching: 

540 # table_source[b] <-> table_cat[use_band[a]] 

541 obj_index = table_source['obj_index'][:] 

542 a, b = esutil.numpy_util.match(use_band, obj_index) 

543 

544 # Update indexes and cut to band-selected stars/sources 

545 table_source['obj_index'][b] = a 

546 _, index_new = np.unique(a, return_index=True) 

547 table_cat[f'source_cat_index_{band}'][use_band] = index_new 

548 

549 # After the following cuts, the catalogs have the following properties: 

550 # - table_cat only contains isolated stars that have at least one source 

551 # in ``band``. 

552 # - table_source only contains ``band`` sources. 

553 # - The slice table_cat["source_cat_index_{band}"]: table_cat["source_cat_index_{band}"] 

554 # + table_cat["nsource_{band}] 

555 # applied to table_source will give all the sources associated with the star. 

556 # - For each source, table_source["obj_index"] points to the index of the associated 

557 # isolated star. 

558 table_source = table_source[b] 

559 table_cat = table_cat[use_band] 

560 

561 # Add reserved flag column to tables 

562 table_cat = np.lib.recfunctions.append_fields( 

563 table_cat, 

564 'reserved', 

565 np.zeros(table_cat.size, dtype=bool), 

566 usemask=False 

567 ) 

568 table_source = np.lib.recfunctions.append_fields( 

569 table_source, 

570 'reserved', 

571 np.zeros(table_source.size, dtype=bool), 

572 usemask=False 

573 ) 

574 

575 # Get reserve star flags 

576 table_cat['reserved'][:] = self.reserve_selection.run( 

577 len(table_cat), 

578 extra=f'{band}_{tract}', 

579 ) 

580 table_source['reserved'][:] = table_cat['reserved'][table_source['obj_index']] 

581 

582 # Offset indexes to account for tract merging 

583 table_cat[f'source_cat_index_{band}'] += merge_source_counter 

584 table_source['obj_index'] += merge_cat_counter 

585 

586 isolated_tables.append(table_cat) 

587 isolated_sources.append(table_source) 

588 

589 merge_cat_counter += len(table_cat) 

590 merge_source_counter += len(table_source) 

591 

592 if len(isolated_tables) > 0: 

593 isolated_table = np.concatenate(isolated_tables) 

594 isolated_source_table = np.concatenate(isolated_sources) 

595 else: 

596 isolated_table = None 

597 isolated_source_table = None 

598 

599 return isolated_table, isolated_source_table 

600 

601 def add_src_colors(self, srcCat, fgcmCat, band): 

602 

603 if self.isPsfDeterminerPiff and fgcmCat is not None: 

604 

605 raSrc = (srcCat['coord_ra'] * u.radian).to(u.degree).value 

606 decSrc = (srcCat['coord_dec'] * u.radian).to(u.degree).value 

607 

608 with Matcher(raSrc, decSrc) as matcher: 

609 idx, idxSrcCat, idxColorCat, d = matcher.query_radius( 

610 fgcmCat["ra"], 

611 fgcmCat["dec"], 

612 1. / 3600.0, 

613 return_indices=True, 

614 ) 

615 

616 magStr1 = self.psf_determiner.config.color[band][0] 

617 magStr2 = self.psf_determiner.config.color[band][2] 

618 colors = fgcmCat[f'mag_{magStr1}'] - fgcmCat[f'mag_{magStr2}'] 

619 

620 for idSrc, idColor in zip(idxSrcCat, idxColorCat): 

621 srcCat[idSrc]['psf_color_value'] = colors[idColor] 

622 srcCat[idSrc]['psf_color_type'] = f"{magStr1}-{magStr2}" 

623 

624 def compute_psf_and_ap_corr_map(self, visit, detector, exposure, src, 

625 isolated_source_table, fgcm_standard_star_cat): 

626 """Compute psf model and aperture correction map for a single exposure. 

627 

628 Parameters 

629 ---------- 

630 visit : `int` 

631 Visit number (for logging). 

632 detector : `int` 

633 Detector number (for logging). 

634 exposure : `lsst.afw.image.ExposureF` 

635 src : `lsst.afw.table.SourceCatalog` 

636 isolated_source_table : `np.ndarray` 

637 

638 Returns 

639 ------- 

640 psf : `lsst.meas.algorithms.ImagePsf` 

641 PSF Model 

642 ap_corr_map : `lsst.afw.image.ApCorrMap` 

643 Aperture correction map. 

644 measured_src : `lsst.afw.table.SourceCatalog` 

645 Updated source catalog with measurements, flags and aperture corrections. 

646 """ 

647 # Extract footprints from the input src catalog for noise replacement. 

648 footprints = SingleFrameMeasurementTask.getFootprintsFromCatalog(src) 

649 

650 # Apply source selector (s/n, flags, etc.) 

651 good_src = self.source_selector.selectSources(src) 

652 if sum(good_src.selected) == 0: 

653 self.log.warning('No good sources remain after cuts for visit %d, detector %d', 

654 visit, detector) 

655 return None, None, None 

656 

657 # Cut down input src to the selected sources 

658 # We use a separate schema/mapper here than for the output/measurement catalog because of 

659 # clashes between fields that were previously run and those that need to be rerun with 

660 # the new psf model. This may be slightly inefficient but keeps input 

661 # and output values cleanly separated. 

662 selection_mapper, selection_schema = self._make_selection_schema_mapper(src.schema) 

663 

664 selected_src = afwTable.SourceCatalog(selection_schema) 

665 selected_src.reserve(good_src.selected.sum()) 

666 selected_src.extend(src[good_src.selected], mapper=selection_mapper) 

667 

668 # The calib flags have been copied from the input table, 

669 # and we reset them here just to ensure they aren't propagated. 

670 selected_src['calib_psf_candidate'] = np.zeros(len(selected_src), dtype=bool) 

671 selected_src['calib_psf_used'] = np.zeros(len(selected_src), dtype=bool) 

672 selected_src['calib_psf_reserved'] = np.zeros(len(selected_src), dtype=bool) 

673 

674 # Find the isolated sources and set flags 

675 matched_src, matched_iso = esutil.numpy_util.match( 

676 selected_src['id'], 

677 isolated_source_table[self.config.id_column] 

678 ) 

679 if len(matched_src) == 0: 

680 self.log.warning( 

681 "No candidates from matched isolate stars for visit=%s, detector=%s " 

682 "(this is probably the result of an earlier astrometry failure).", 

683 visit, detector, 

684 ) 

685 return None, None, None 

686 

687 matched_arr = np.zeros(len(selected_src), dtype=bool) 

688 matched_arr[matched_src] = True 

689 selected_src['calib_psf_candidate'] = matched_arr 

690 

691 reserved_arr = np.zeros(len(selected_src), dtype=bool) 

692 reserved_arr[matched_src] = isolated_source_table['reserved'][matched_iso] 

693 selected_src['calib_psf_reserved'] = reserved_arr 

694 

695 selected_src = selected_src[selected_src['calib_psf_candidate']].copy(deep=True) 

696 

697 # Make the measured source catalog as well, based on the selected catalog. 

698 measured_src = afwTable.SourceCatalog(self.schema) 

699 measured_src.reserve(len(selected_src)) 

700 measured_src.extend(selected_src, mapper=self.schema_mapper) 

701 

702 # We need to copy over the calib_psf flags because they were not in the mapper 

703 measured_src['calib_psf_candidate'] = selected_src['calib_psf_candidate'] 

704 measured_src['calib_psf_reserved'] = selected_src['calib_psf_reserved'] 

705 if exposure.filter.hasBandLabel(): 

706 band = exposure.filter.bandLabel 

707 else: 

708 band = None 

709 self.add_src_colors(selected_src, fgcm_standard_star_cat, band) 

710 self.add_src_colors(measured_src, fgcm_standard_star_cat, band) 

711 

712 # Select the psf candidates from the selection catalog 

713 try: 

714 psf_selection_result = self.make_psf_candidates.run(selected_src, exposure=exposure) 

715 _ = self.make_psf_candidates.run(measured_src, exposure=exposure) 

716 except Exception as e: 

717 self.log.exception('Failed to make PSF candidates for visit %d, detector %d: %s', 

718 visit, detector, e) 

719 return None, None, measured_src 

720 

721 psf_cand_cat = psf_selection_result.goodStarCat 

722 

723 # Make list of psf candidates to send to the determiner 

724 # (omitting those marked as reserved) 

725 psf_determiner_list = [cand for cand, use 

726 in zip(psf_selection_result.psfCandidates, 

727 ~psf_cand_cat['calib_psf_reserved']) if use] 

728 flag_key = psf_cand_cat.schema['calib_psf_used'].asKey() 

729 

730 try: 

731 psf, cell_set = self.psf_determiner.determinePsf(exposure, 

732 psf_determiner_list, 

733 self.metadata, 

734 flagKey=flag_key) 

735 except Exception as e: 

736 self.log.exception('Failed to determine PSF for visit %d, detector %d: %s', 

737 visit, detector, e) 

738 return None, None, measured_src 

739 # Verify that the PSF is usable by downstream tasks 

740 sigma = psf.computeShape(psf.getAveragePosition(), psf.getAverageColor()).getDeterminantRadius() 

741 if np.isnan(sigma): 

742 self.log.warning('Failed to determine psf for visit %d, detector %d: ' 

743 'Computed final PSF size is NAN.', 

744 visit, detector) 

745 return None, None, measured_src 

746 

747 # Set the psf in the exposure for measurement/aperture corrections. 

748 exposure.setPsf(psf) 

749 

750 # At this point, we need to transfer the psf used flag from the selection 

751 # catalog to the measurement catalog. 

752 matched_selected, matched_measured = esutil.numpy_util.match( 

753 selected_src['id'], 

754 measured_src['id'] 

755 ) 

756 measured_used = np.zeros(len(measured_src), dtype=bool) 

757 measured_used[matched_measured] = selected_src['calib_psf_used'][matched_selected] 

758 measured_src['calib_psf_used'] = measured_used 

759 

760 # Next, we do the measurement on all the psf candidate, used, and reserved stars. 

761 # We use the full footprint list from the input src catalog for noise replacement. 

762 try: 

763 self.measurement.run(measCat=measured_src, exposure=exposure, footprints=footprints) 

764 except Exception as e: 

765 self.log.warning('Failed to make measurements for visit %d, detector %d: %s', 

766 visit, detector, e) 

767 return psf, None, measured_src 

768 

769 # And finally the ap corr map. 

770 try: 

771 ap_corr_map = self.measure_ap_corr.run(exposure=exposure, 

772 catalog=measured_src).apCorrMap 

773 except Exception as e: 

774 self.log.warning('Failed to compute aperture corrections for visit %d, detector %d: %s', 

775 visit, detector, e) 

776 return psf, None, measured_src 

777 

778 # Need to merge the original normalization aperture correction map. 

779 ap_corr_map_input = exposure.apCorrMap 

780 for key in ap_corr_map_input: 

781 if key not in ap_corr_map: 

782 ap_corr_map[key] = ap_corr_map_input[key] 

783 

784 self.apply_ap_corr.run(catalog=measured_src, apCorrMap=ap_corr_map) 

785 

786 return psf, ap_corr_map, measured_src 

787 

788 

789class FinalizeCharacterizationTask(FinalizeCharacterizationTaskBase): 

790 """Run final characterization on full visits.""" 

791 ConfigClass = FinalizeCharacterizationConfig 

792 _DefaultName = 'finalize_characterization' 

793 

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

795 input_handle_dict = butlerQC.get(inputRefs) 

796 

797 band = butlerQC.quantum.dataId['band'] 

798 visit = butlerQC.quantum.dataId['visit'] 

799 

800 src_dict_temp = {handle.dataId['detector']: handle 

801 for handle in input_handle_dict['srcs']} 

802 calexp_dict_temp = {handle.dataId['detector']: handle 

803 for handle in input_handle_dict['calexps']} 

804 isolated_star_cat_dict_temp = {handle.dataId['tract']: handle 

805 for handle in input_handle_dict['isolated_star_cats']} 

806 isolated_star_source_dict_temp = {handle.dataId['tract']: handle 

807 for handle in input_handle_dict['isolated_star_sources']} 

808 

809 src_dict = {detector: src_dict_temp[detector] for 

810 detector in sorted(src_dict_temp.keys())} 

811 calexp_dict = {detector: calexp_dict_temp[detector] for 

812 detector in sorted(calexp_dict_temp.keys())} 

813 isolated_star_cat_dict = {tract: isolated_star_cat_dict_temp[tract] for 

814 tract in sorted(isolated_star_cat_dict_temp.keys())} 

815 isolated_star_source_dict = {tract: isolated_star_source_dict_temp[tract] for 

816 tract in sorted(isolated_star_source_dict_temp.keys())} 

817 

818 if self.config.psf_determiner['piff'].useColor: 

819 fgcm_standard_star_dict_temp = {handle.dataId['tract']: handle 

820 for handle in input_handle_dict['fgcm_standard_star']} 

821 fgcm_standard_star_dict = {tract: fgcm_standard_star_dict_temp[tract] for 

822 tract in sorted(fgcm_standard_star_dict_temp.keys())} 

823 else: 

824 fgcm_standard_star_dict = None 

825 

826 struct = self.run( 

827 visit, 

828 band, 

829 isolated_star_cat_dict, 

830 isolated_star_source_dict, 

831 src_dict, 

832 calexp_dict, 

833 fgcm_standard_star_dict=fgcm_standard_star_dict, 

834 ) 

835 

836 butlerQC.put(struct.psf_ap_corr_cat, outputRefs.finalized_psf_ap_corr_cat) 

837 butlerQC.put(struct.output_table, outputRefs.finalized_src_table) 

838 

839 def run(self, visit, band, isolated_star_cat_dict, isolated_star_source_dict, 

840 src_dict, calexp_dict, fgcm_standard_star_dict=None): 

841 """ 

842 Run the FinalizeCharacterizationTask. 

843 

844 Parameters 

845 ---------- 

846 visit : `int` 

847 Visit number. Used in the output catalogs. 

848 band : `str` 

849 Band name. Used to select reserved stars. 

850 isolated_star_cat_dict : `dict` 

851 Per-tract dict of isolated star catalog handles.