Coverage for python / lsst / fgcmcal / fgcmBuildFromIsolatedStars.py: 0%

1# See COPYRIGHT file at the top of the source tree. 

2# 

3# This file is part of fgcmcal. 

4# 

5# Developed for the LSST Data Management System. 

6# This product includes software developed by the LSST Project 

7# (https://www.lsst.org). 

8# See the COPYRIGHT file at the top-level directory of this distribution 

9# for details of code ownership. 

10# 

11# This program is free software: you can redistribute it and/or modify 

12# it under the terms of the GNU General Public License as published by 

13# the Free Software Foundation, either version 3 of the License, or 

14# (at your option) any later version. 

15# 

16# This program is distributed in the hope that it will be useful, 

17# but WITHOUT ANY WARRANTY; without even the implied warranty of 

18# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

19# GNU General Public License for more details. 

20# 

21# You should have received a copy of the GNU General Public License 

22# along with this program. If not, see <https://www.gnu.org/licenses/>. 

23"""Build star observations for input to FGCM using sourceTable_visit. 

24 

25This task finds all the visits and sourceTable_visits in a repository (or a 

26subset based on command line parameters) and extracts all the potential 

27calibration stars for input into fgcm. This task additionally uses fgcm to 

28match star observations into unique stars, and performs as much cleaning of the 

29input catalog as possible. 

30""" 

31import warnings 

32import numpy as np 

33import esutil 

34import hpgeom as hpg 

35from smatch.matcher import Matcher 

36from astropy.table import Table, vstack 

37 

38from fgcm.fgcmUtilities import objFlagDict, histogram_rev_sorted, getMemoryString 

39 

40import lsst.afw.cameraGeom as afwCameraGeom 

41import lsst.pex.config as pexConfig 

42import lsst.pipe.base as pipeBase 

43from lsst.pipe.base import connectionTypes 

44from lsst.meas.algorithms import ReferenceObjectLoader, LoadReferenceObjectsConfig 

45from lsst.pipe.tasks.reserveIsolatedStars import ReserveIsolatedStarsTask 

46 

47from .fgcmBuildStarsBase import FgcmBuildStarsConfigBase, FgcmBuildStarsBaseTask 

48from .utilities import computeApproxPixelAreaFields, computeApertureRadiusFromName 

49 

50__all__ = ["FgcmBuildFromIsolatedStarsConfig", "FgcmBuildFromIsolatedStarsTask"] 

51 

52 

53class FgcmBuildFromIsolatedStarsConnections(pipeBase.PipelineTaskConnections, 

54 dimensions=("instrument",), 

55 defaultTemplates={}): 

56 camera = connectionTypes.PrerequisiteInput( 

57 doc="Camera instrument", 

58 name="camera", 

59 storageClass="Camera", 

60 dimensions=("instrument",), 

61 isCalibration=True, 

62 ) 

63 fgcm_lookup_table = connectionTypes.PrerequisiteInput( 

64 doc=("Atmosphere + instrument look-up-table for FGCM throughput and " 

65 "chromatic corrections."), 

66 name="fgcmLookUpTable", 

67 storageClass="Catalog", 

68 dimensions=("instrument",), 

69 deferLoad=True, 

70 ) 

71 ref_cat = connectionTypes.PrerequisiteInput( 

72 doc="Reference catalog to use for photometric calibration.", 

73 name="cal_ref_cat", 

74 storageClass="SimpleCatalog", 

75 dimensions=("skypix",), 

76 deferLoad=True, 

77 multiple=True, 

78 ) 

79 isolated_star_cats = pipeBase.connectionTypes.Input( 

80 doc=("Catalog of isolated stars with average positions, number of associated " 

81 "sources, and indexes to the isolated_star_sources catalogs."), 

82 name="isolated_star_presource_associations", 

83 storageClass="ArrowAstropy", 

84 dimensions=("instrument", "tract", "skymap"), 

85 deferLoad=True, 

86 multiple=True, 

87 ) 

88 isolated_star_sources = pipeBase.connectionTypes.Input( 

89 doc=("Catalog of isolated star sources with sourceIds, and indexes to the " 

90 "isolated_star_cats catalogs."), 

91 name="isolated_star_presources", 

92 storageClass="ArrowAstropy", 

93 dimensions=("instrument", "tract", "skymap"), 

94 deferLoad=True, 

95 multiple=True, 

96 ) 

97 visit_summaries = connectionTypes.Input( 

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

99 "detector id for the id and must be sorted for fast lookups of a " 

100 "detector."), 

101 name="visitSummary", 

102 storageClass="ExposureCatalog", 

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

104 deferLoad=True, 

105 multiple=True, 

106 ) 

107 fgcm_visit_catalog = connectionTypes.Output( 

108 doc="Catalog of visit information for fgcm.", 

109 name="fgcmVisitCatalog", 

110 storageClass="Catalog", 

111 dimensions=("instrument",), 

112 ) 

113 fgcm_star_observations = connectionTypes.Output( 

114 doc="Catalog of star observations for fgcm.", 

115 name="fgcm_star_observations", 

116 storageClass="ArrowAstropy", 

117 dimensions=("instrument",), 

118 ) 

119 fgcm_star_ids = connectionTypes.Output( 

120 doc="Catalog of fgcm calibration star IDs.", 

121 name="fgcm_star_ids", 

122 storageClass="ArrowAstropy", 

123 dimensions=("instrument",), 

124 ) 

125 fgcm_reference_stars = connectionTypes.Output( 

126 doc="Catalog of fgcm-matched reference stars.", 

127 name="fgcm_reference_stars", 

128 storageClass="ArrowAstropy", 

129 dimensions=("instrument",), 

130 ) 

131 

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

133 super().__init__(config=config) 

134 

135 if not config.doReferenceMatches: 

136 self.prerequisiteInputs.remove("ref_cat") 

137 self.outputs.remove("fgcm_reference_stars") 

138 

139 def getSpatialBoundsConnections(self): 

140 return ("isolated_star_cats", "visit_summaries") 

141 

142 

143class FgcmBuildFromIsolatedStarsConfig(FgcmBuildStarsConfigBase, pipeBase.PipelineTaskConfig, 

144 pipelineConnections=FgcmBuildFromIsolatedStarsConnections): 

145 """Config for FgcmBuildFromIsolatedStarsTask.""" 

146 referenceCCD = pexConfig.Field( 

147 doc="Reference detector for checking PSF and background.", 

148 dtype=int, 

149 default=40, 

150 ) 

151 reserve_selection = pexConfig.ConfigurableField( 

152 target=ReserveIsolatedStarsTask, 

153 doc="Task to select reserved stars.", 

154 ) 

155 

156 def setDefaults(self): 

157 super().setDefaults() 

158 

159 self.reserve_selection.reserve_name = "fgcmcal" 

160 self.reserve_selection.reserve_fraction = 0.1 

161 

162 # The names here correspond to the isolated_star_sources. 

163 self.instFluxField = 'normCompTophatFlux_instFlux' 

164 self.localBackgroundFluxField = 'localBackground_instFlux' 

165 self.apertureInnerInstFluxField = 'apFlux_12_0_instFlux' 

166 self.apertureOuterInstFluxField = 'apFlux_17_0_instFlux' 

167 

168 source_selector = self.sourceSelector["science"] 

169 source_selector.setDefaults() 

170 

171 source_selector.doFlags = False 

172 source_selector.doSignalToNoise = True 

173 source_selector.doUnresolved = False 

174 source_selector.doIsolated = False 

175 source_selector.doRequireFiniteRaDec = False 

176 

177 source_selector.flags.bad = [] 

178 

179 source_selector.signalToNoise.minimum = 9.0 

180 source_selector.signalToNoise.maximum = 1000.0 

181 source_selector.signalToNoise.fluxField = self.instFluxField 

182 source_selector.signalToNoise.errField = self.instFluxField + 'Err' 

183 

184 

185class FgcmBuildFromIsolatedStarsTask(FgcmBuildStarsBaseTask): 

186 """Build star catalog for FGCM global calibration, using the isolated star catalogs. 

187 """ 

188 ConfigClass = FgcmBuildFromIsolatedStarsConfig 

189 _DefaultName = "fgcmBuildFromIsolatedStars" 

190 

191 canMultiprocess = False 

192 

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

194 super().__init__(**kwargs) 

195 self.makeSubtask('reserve_selection') 

196 

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

198 input_ref_dict = butlerQC.get(inputRefs) 

199 

200 isolated_star_cat_handles = input_ref_dict["isolated_star_cats"] 

201 isolated_star_source_handles = input_ref_dict["isolated_star_sources"] 

202 

203 isolated_star_cat_handle_dict = { 

204 handle.dataId["tract"]: handle for handle in isolated_star_cat_handles 

205 } 

206 isolated_star_source_handle_dict = { 

207 handle.dataId["tract"]: handle for handle in isolated_star_source_handles 

208 } 

209 

210 if len(isolated_star_cat_handle_dict) != len(isolated_star_source_handle_dict): 

211 raise RuntimeError("isolated_star_cats and isolate_star_sources must have same length.") 

212 

213 for tract in isolated_star_cat_handle_dict: 

214 if tract not in isolated_star_source_handle_dict: 

215 raise RuntimeError(f"tract {tract} in isolated_star_cats but not isolated_star_sources") 

216 

217 if self.config.doReferenceMatches: 

218 

219 # Prepare the reference catalog loader 

220 ref_config = LoadReferenceObjectsConfig() 

221 ref_config.filterMap = self.config.fgcmLoadReferenceCatalog.filterMap 

222 ref_obj_loader = ReferenceObjectLoader(dataIds=[ref.datasetRef.dataId 

223 for ref in inputRefs.ref_cat], 

224 refCats=butlerQC.get(inputRefs.ref_cat), 

225 name=self.config.connections.ref_cat, 

226 log=self.log, 

227 config=ref_config) 

228 self.makeSubtask('fgcmLoadReferenceCatalog', 

229 refObjLoader=ref_obj_loader, 

230 refCatName=self.config.connections.ref_cat) 

231 

232 # The visit summary handles for use with fgcmMakeVisitCatalog must be keyed with 

233 # visit, and values are a list with the first value being the visit_summary_handle, 

234 # the second is the source catalog (which is not used, but a place-holder is 

235 # left for compatibility reasons.) 

236 visit_summary_handle_dict = {handle.dataId['visit']: [handle, None] for 

237 handle in input_ref_dict['visit_summaries']} 

238 

239 camera = input_ref_dict["camera"] 

240 

241 lookup_table_handle = input_ref_dict["fgcm_lookup_table"] 

242 

243 struct = self.run( 

244 camera=camera, 

245 visit_summary_handle_dict=visit_summary_handle_dict, 

246 isolated_star_cat_handle_dict=isolated_star_cat_handle_dict, 

247 isolated_star_source_handle_dict=isolated_star_source_handle_dict, 

248 lookup_table_handle=lookup_table_handle, 

249 ) 

250 

251 butlerQC.put(struct.fgcm_visit_catalog, outputRefs.fgcm_visit_catalog) 

252 butlerQC.put(struct.fgcm_star_observations, outputRefs.fgcm_star_observations) 

253 butlerQC.put(struct.fgcm_star_ids, outputRefs.fgcm_star_ids) 

254 if self.config.doReferenceMatches: 

255 butlerQC.put(struct.fgcm_reference_stars, outputRefs.fgcm_reference_stars) 

256 

257 def run(self, *, camera, visit_summary_handle_dict, isolated_star_cat_handle_dict, 

258 isolated_star_source_handle_dict, lookup_table_handle): 

259 """Run the fgcmBuildFromIsolatedStarsTask. 

260 

261 Parameters 

262 ---------- 

263 camera : `lsst.afw.cameraGeom.Camera` 

264 Camera object. 

265 visit_summary_handle_dict : `dict` [`int`, [`lsst.daf.butler.DeferredDatasetHandle`]] 

266 Visit summary dataset handles, with the visit as key. 

267 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`] 

268 Isolated star catalog dataset handles, with the tract as key. 

269 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`] 

270 Isolated star source dataset handles, with the tract as key. 

271 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle` 

272 Data reference to fgcm look-up table. 

273 

274 Returns 

275 ------- 

276 struct : `lsst.pipe.base.struct` 

277 Catalogs for persistence, with attributes: 

278 

279 ``fgcm_visit_catalog`` 

280 Catalog of per-visit data (`lsst.afw.table.ExposureCatalog`). 

281 ``fgcm_star_observations`` 

282 Catalog of individual star observations (`astropy.table.Table`). 

283 ``fgcm_star_ids`` 

284 Catalog of per-star ids and positions (`astropy.table.Table`). 

285 ``fgcm_reference_stars`` 

286 Catalog of reference stars matched to fgcm stars (`astropy.table.Table`). 

287 """ 

288 # Compute aperture radius if necessary. This is useful to do now before 

289 # any heave lifting has happened (fail early). 

290 calib_flux_aperture_radius = None 

291 if self.config.doSubtractLocalBackground: 

292 try: 

293 calib_flux_aperture_radius = computeApertureRadiusFromName(self.config.instFluxField) 

294 except RuntimeError as e: 

295 raise RuntimeError("Could not determine aperture radius from %s. " 

296 "Cannot use doSubtractLocalBackground." % 

297 (self.config.instFluxField)) from e 

298 

299 # Check that we have the lookup_table_handle if we are doing reference matches. 

300 if self.config.doReferenceMatches: 

301 if lookup_table_handle is None: 

302 raise RuntimeError("Must supply lookup_table_handle if config.doReferenceMatches is True.") 

303 

304 lut_cat = lookup_table_handle.get() 

305 if len(camera) == lut_cat[0]["nCcd"]: 

306 use_science_detectors = False 

307 else: 

308 # If the LUT has a different number of detectors than 

309 # the camera, then we only want to use science detectors 

310 # in the focal plane projector. 

311 # Note that in the LUT building code we either have 

312 # all detectors or only science detectors. 

313 use_science_detectors = True 

314 del lut_cat 

315 

316 visit_cat = self.fgcmMakeVisitCatalog( 

317 camera, 

318 visit_summary_handle_dict, 

319 useScienceDetectors=use_science_detectors, 

320 ) 

321 

322 fgcm_obj = self._make_all_star_obs_from_isolated_stars( 

323 isolated_star_cat_handle_dict, 

324 isolated_star_source_handle_dict, 

325 visit_cat, 

326 camera, 

327 calib_flux_aperture_radius=calib_flux_aperture_radius, 

328 use_science_detectors=use_science_detectors, 

329 ) 

330 

331 star_ids = self._density_downsample(fgcm_obj) 

332 

333 # Select and concatenate the isolated stars and sources. 

334 fgcm_obj, star_obs = self._make_all_star_obs_from_isolated_stars( 

335 isolated_star_cat_handle_dict, 

336 isolated_star_source_handle_dict, 

337 visit_cat, 

338 camera, 

339 calib_flux_aperture_radius=calib_flux_aperture_radius, 

340 use_science_detectors=use_science_detectors, 

341 input_star_ids=star_ids, 

342 return_observations=True, 

343 ) 

344 

345 self._compute_delta_aper_summary_statistics( 

346 visit_cat, 

347 star_obs, 

348 ) 

349 

350 # Mark reserve stars 

351 self._mark_reserve_stars(fgcm_obj) 

352 

353 # Do reference association. 

354 if self.config.doReferenceMatches: 

355 fgcm_ref = self._associate_reference_stars(lookup_table_handle, visit_cat, fgcm_obj) 

356 else: 

357 fgcm_ref = None 

358 

359 return pipeBase.Struct( 

360 fgcm_visit_catalog=visit_cat, 

361 fgcm_star_observations=star_obs, 

362 fgcm_star_ids=fgcm_obj, 

363 fgcm_reference_stars=fgcm_ref, 

364 ) 

365 

366 def _make_all_star_obs_from_isolated_stars( 

367 self, 

368 isolated_star_cat_handle_dict, 

369 isolated_star_source_handle_dict, 

370 visit_cat, 

371 camera, 

372 calib_flux_aperture_radius=None, 

373 use_science_detectors=False, 

374 input_star_ids=None, 

375 return_observations=False, 

376 ): 

377 """ 

378 Make all the star observations from isolated star catalogs. 

379 

380 This is intended to be run in two passes. In the first pass, a limited 

381 set of columns is read in (only what is required to select potential 

382 "good stars") and the stars are returned. In the second pass, the full 

383 set of columns is read and only stars that are in input_star_ids 

384 are put into the aggregated tables. 

385 

386 Parameters 

387 ---------- 

388 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`] 

389 Isolated star catalog dataset handles, with the tract as key. 

390 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`] 

391 Isolated star source dataset handles, with the tract as key. 

392 visit_cat : `lsst.afw.table.ExposureCatalog` 

393 Catalog of per-visit data. 

394 camera : `lsst.afw.cameraGeom.Camera` 

395 The camera object. 

396 calib_flux_aperture_radius : `float`, optional 

397 Radius for the calibration flux aperture. 

398 use_science_detectors : `bool`, optional 

399 Only use science detectors? 

400 input_star_ids : `np.ndarray`, optional 

401 Array of star IDs to use. Used for second pass and 

402 full load. 

403 return_observations : `bool`, optional 

404 Return the observations? Should be False for first pass. 

405 

406 Returns 

407 ------- 

408 fgcm_obj : `astropy.table.Table` 

409 Catalog of ids and positions for each star. 

410 star_obs : `astropy.table.Table` 

411 Catalog of individual star observations if 

412 return_observations is True. 

413 """ 

414 # Ensure that input stars are sorted. 

415 if input_star_ids is not None: 

416 input_star_ids = np.sort(input_star_ids) 

417 

418 source_columns = [ 

419 "sourceId", 

420 "visit", 

421 "detector", 

422 "physical_filter", 

423 "band", 

424 "obj_index", 

425 self.config.instFluxField, 

426 self.config.instFluxField + "Err", 

427 ] 

428 

429 if return_observations: 

430 # We only need these columns if we are returning all 

431 # the observations. 

432 source_columns.extend( 

433 [ 

434 "ra", 

435 "dec", 

436 "x", 

437 "y", 

438 self.config.apertureInnerInstFluxField, 

439 self.config.apertureInnerInstFluxField + "Err", 

440 self.config.apertureOuterInstFluxField, 

441 self.config.apertureOuterInstFluxField + "Err", 

442 ], 

443 ) 

444 

445 if self.config.doSubtractLocalBackground: 

446 source_columns.append(self.config.localBackgroundFluxField) 

447 local_background_flag_name = self.config.localBackgroundFluxField[0: -len('instFlux')] + 'flag' 

448 source_columns.append(local_background_flag_name) 

449 

450 if self.sourceSelector.config.doFlags: 

451 source_columns.extend(self.sourceSelector.config.flags.bad) 

452 

453 if self.config.doSubtractLocalBackground: 

454 local_background_area = np.pi*calib_flux_aperture_radius**2. 

455 

456 # Compute the approximate pixel area over the full focal plane 

457 # from the WCS jacobian using the camera model. 

458 approx_pixel_area_fields = computeApproxPixelAreaFields(camera) 

459 

460 # Construct mapping from detector number to index. 

461 detector_mapping = {} 

462 for detector_index, detector in enumerate(camera): 

463 detector_mapping[detector.getId()] = detector_index 

464 

465 fgcm_objs = [] 

466 star_obs_cats = [] 

467 merge_source_counter = 0 

468 

469 k = 2.5/np.log(10.) 

470 

471 visit_cat_table = visit_cat.asAstropy() 

472 

473 n_tracts = len(isolated_star_cat_handle_dict) 

474 

475 for index, tract in enumerate(sorted(isolated_star_cat_handle_dict)): 

476 if (index % 50) == 0: 

477 self.log.info(getMemoryString(f"Loading tract index {index}")) 

478 

479 stars = isolated_star_cat_handle_dict[tract].get() 

480 sources = isolated_star_source_handle_dict[tract].get(parameters={"columns": source_columns}) 

481 

482 # Down-select sources. 

483 good_sources = self.sourceSelector.selectSources(sources).selected 

484 if self.config.doSubtractLocalBackground: 

485 good_sources &= (~sources[local_background_flag_name]) 

486 local_background = local_background_area*sources[self.config.localBackgroundFluxField] 

487 good_sources &= ((sources[self.config.instFluxField] - local_background) > 0) 

488 

489 # We need to do this check as early as possible to 

490 # ensure we can do the match test below. 

491 if good_sources.sum() == 0: 

492 self.log.info("No good sources found in tract %d", tract) 

493 continue 

494 

495 # We also need to make sure that all sources are in the 

496 # input list of visits and detectors. 

497 

498 detector_ids = [] 

499 for detector in camera: 

500 if use_science_detectors: 

501 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE: 

502 continue 

503 detector_ids.append(detector.getId()) 

504 

505 missing_sources1 = np.ones(len(sources), dtype=np.bool_) 

506 _, sources_match = esutil.numpy_util.match(visit_cat_table["visit"], sources["visit"]) 

507 missing_sources1[sources_match] = False 

508 good_sources &= (~missing_sources1) 

509 

510 missing_sources2 = np.ones(len(sources), dtype=np.bool_) 

511 _, sources_match = esutil.numpy_util.match(detector_ids, sources["detector"]) 

512 missing_sources2[sources_match] = False 

513 good_sources &= (~missing_sources2) 

514 

515 # Check again after further tests. 

516 if good_sources.sum() == 0: 

517 self.log.info("No good sources found in tract %d", tract) 

518 continue 

519 

520 self.log.info( 

521 "Tract %d (%d/%d) contains %d good sources.", 

522 tract, 

523 index, 

524 n_tracts, 

525 good_sources.sum(), 

526 ) 

527 

528 # Need to count the observations of each star after cuts, per band. 

529 # If we have requiredBands specified, we must make sure that each star 

530 # has the minumum number of observations in each of thos bands. 

531 # Otherwise, we must make sure that each star has at least the minimum 

532 # number of observations in _any_ band. 

533 if len(self.config.requiredBands) > 0: 

534 loop_bands = self.config.requiredBands 

535 else: 

536 loop_bands = np.unique(sources["band"]) 

537 

538 n_req = np.zeros((len(loop_bands), len(stars)), dtype=np.int32) 

539 for i, band in enumerate(loop_bands): 

540 (band_use,) = (sources[good_sources]["band"] == band).nonzero() 

541 np.add.at( 

542 n_req, 

543 (i, sources[good_sources]["obj_index"][band_use]), 

544 1, 

545 ) 

546 

547 if len(self.config.requiredBands) > 0: 

548 # The min gives us the band with the fewest observations, which must be 

549 # above the limit. 

550 (good_stars,) = (n_req.min(axis=0) >= self.config.minPerBand).nonzero() 

551 else: 

552 # The max gives us the band with the most observations, which must be 

553 # above the limit. 

554 (good_stars,) = (n_req.max(axis=0) >= self.config.minPerBand).nonzero() 

555 

556 if len(good_stars) == 0: 

557 self.log.info("No good stars found in tract %d", tract) 

558 continue 

559 

560 if not return_observations: 

561 fgcm_obj = Table() 

562 fgcm_obj["isolated_star_id"] = stars["isolated_star_id"][good_stars] 

563 fgcm_obj["ra"] = stars["ra"][good_stars] 

564 fgcm_obj["dec"] = stars["dec"][good_stars] 

565 fgcm_obj["fgcm_id"] = np.zeros(len(fgcm_obj), dtype=np.int64) 

566 

567 fgcm_objs.append(fgcm_obj) 

568 

569 continue 

570 elif input_star_ids is not None: 

571 # Down-select good_stars indices to those in the inputs. 

572 sub1 = np.clip( 

573 np.searchsorted(input_star_ids, stars["isolated_star_id"][good_stars]), 

574 0, 

575 len(input_star_ids) - 1, 

576 ) 

577 matched = (input_star_ids[sub1] == stars["isolated_star_id"][good_stars]) 

578 good_stars = good_stars[matched] 

579 

580 if len(good_stars) == 0: 

581 self.log.info("No good stars found after down-selection in tract %d", tract) 

582 continue 

583 

584 # The following is only used if we are assembling the observations. 

585 

586 # With the following matching: 

587 # sources[good_sources][b] <-> stars[good_stars[a]] 

588 obj_index = sources["obj_index"][good_sources] 

589 a, b = esutil.numpy_util.match(good_stars, obj_index) 

590 

591 # Update indexes and cut to selected stars/sources 

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

593 stars["source_cat_index"][good_stars] = index_new 

594 sources = sources[good_sources][b] 

595 sources["obj_index"][:] = a 

596 stars = stars[good_stars] 

597 

598 nsource = np.zeros(len(stars), dtype=np.int32) 

599 np.add.at( 

600 nsource, 

601 sources["obj_index"], 

602 1, 

603 ) 

604 stars["nsource"][:] = nsource 

605 

606 # After these cuts, the catalogs have the following properties: 

607 # - ``stars`` only contains isolated stars that have the minimum number of good 

608 # sources in the required bands. 

609 # - ``sources`` has been cut to the good sources. 

610 # - The slice [stars["source_cat_index"]: stars["source_cat_index"] 

611 # + stars["nsource"]] 

612 # applied to ``sources`` will give all the sources associated with the star. 

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

614 # isolated star. 

615 

616 # We now reformat the sources and compute the ``observations`` that fgcm expects. 

617 star_obs = Table() 

618 star_obs["ra"] = sources["ra"].astype(np.float64) 

619 star_obs["dec"] = sources["dec"].astype(np.float64) 

620 star_obs["x"] = sources["x"].astype(np.float64) 

621 star_obs["y"] = sources["y"].astype(np.float64) 

622 star_obs["visit"] = sources["visit"].astype(np.int64) 

623 star_obs["detector"] = sources["detector"].astype(np.int32) 

624 star_obs["inst_mag"] = np.zeros(len(sources), dtype="f4") 

625 star_obs["inst_mag_err"] = np.zeros(len(sources), dtype="f4") 

626 star_obs["jacobian"] = np.zeros(len(sources), dtype="f4") 

627 star_obs["delta_mag_bkg"] = np.zeros(len(sources), dtype="f4") 

628 star_obs["delta_mag_aper"] = np.zeros(len(sources), dtype="f4") 

629 star_obs["delta_mag_err_aper"] = np.zeros(len(sources), dtype="f4") 

630 

631 visit_match, obs_match = esutil.numpy_util.match(visit_cat_table["visit"], sources["visit"]) 

632 

633 exp_time = np.zeros(len(star_obs)) 

634 exp_time[obs_match] = visit_cat_table["exptime"][visit_match] 

635 

636 with warnings.catch_warnings(): 

637 # Ignore warnings, we will filter infinities and nans below. 

638 warnings.simplefilter("ignore") 

639 

640 inst_mag_inner = -2.5*np.log10(sources[self.config.apertureInnerInstFluxField]) 

641 inst_mag_err_inner = k*(sources[self.config.apertureInnerInstFluxField + "Err"] 

642 / sources[self.config.apertureInnerInstFluxField]) 

643 inst_mag_outer = -2.5*np.log10(sources[self.config.apertureOuterInstFluxField]) 

644 inst_mag_err_outer = k*(sources[self.config.apertureOuterInstFluxField + "Err"] 

645 / sources[self.config.apertureOuterInstFluxField]) 

646 star_obs["delta_mag_aper"][:] = inst_mag_inner - inst_mag_outer 

647 star_obs["delta_mag_err_aper"][:] = np.sqrt(inst_mag_err_inner**2. + inst_mag_err_outer**2.) 

648 # Set bad values to sentinel values for fgcm. 

649 bad = ~np.isfinite(star_obs["delta_mag_aper"]) 

650 star_obs["delta_mag_aper"][bad] = 99.0 

651 star_obs["delta_mag_err_aper"][bad] = 99.0 

652 

653 if self.config.doSubtractLocalBackground: 

654 # At the moment we only adjust the flux and not the flux 

655 # error by the background because the error on 

656 # base_LocalBackground_instFlux is the rms error in the 

657 # background annulus, not the error on the mean in the 

658 # background estimate (which is much smaller, by sqrt(n) 

659 # pixels used to estimate the background, which we do not 

660 # have access to in this task). In the default settings, 

661 # the annulus is sufficiently large such that these 

662 # additional errors are negligibly small (much less 

663 # than a mmag in quadrature). 

664 

665 # This is the difference between the mag with local background correction 

666 # and the mag without local background correction. 

667 local_background = local_background_area*sources[self.config.localBackgroundFluxField] 

668 star_obs["delta_mag_bkg"][:] = ( 

669 -2.5*np.log10(sources[self.config.instFluxField] 

670 - local_background) - 

671 -2.5*np.log10(sources[self.config.instFluxField]) 

672 ) 

673 

674 # Need to loop over detectors here. 

675 for detector in camera: 

676 if use_science_detectors: 

677 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE: 

678 continue 

679 

680 detector_id = detector.getId() 

681 # used index for all observations with a given detector 

682 (use,) = (star_obs["detector"][obs_match] == detector_id).nonzero() 

683 # Prior to running the calibration, we want to remove the effect 

684 # of the jacobian of the WCS because that is a known quantity. 

685 # Ideally, this would be done for each individual WCS, but this 

686 # is extremely slow and makes small differences that are much 

687 # smaller than the variation in the throughput due to other issues. 

688 # Therefore, we use the approximate jacobian estimated from the 

689 # camera model. 

690 jac = approx_pixel_area_fields[detector_id].evaluate( 

691 star_obs["x"][obs_match][use], 

692 star_obs["y"][obs_match][use], 

693 ) 

694 star_obs["jacobian"][obs_match[use]] = jac 

695 scaled_inst_flux = (sources[self.config.instFluxField][obs_match[use]] 

696 * visit_cat_table["scaling"][visit_match[use], 

697 detector_mapping[detector_id]]) 

698 star_obs["inst_mag"][obs_match[use]] = (-2.5 * np.log10(scaled_inst_flux 

699 / exp_time[use])) 

700 

701 # Compute instMagErr from inst_flux_err/inst_flux. 

702 # This does not need to be computed per-detector because there 

703 # is no need for the per-detector scaling which cancels out. 

704 star_obs["inst_mag_err"][:] = k*(sources[self.config.instFluxField + "Err"] 

705 / sources[self.config.instFluxField]) 

706 

707 # Apply the jacobian if configured to do so. 

708 if self.config.doApplyWcsJacobian: 

709 star_obs["inst_mag"][:] -= 2.5*np.log10(star_obs["jacobian"]) 

710 

711 # We now reformat the stars and compute the ''objects'' that fgcm expects. 

712 # Note that we are very careful here about casting to the correct 

713 # data types. In particular, obs_arr_index must be a 64 bit integer 

714 # even if each individual isolated_star_sources can use a 32 bit 

715 # index counter. 

716 fgcm_obj = Table() 

717 fgcm_obj["fgcm_id"] = np.zeros(len(stars), dtype=np.int64) 

718 fgcm_obj["isolated_star_id"] = stars["isolated_star_id"].astype(np.int64) 

719 fgcm_obj["ra"] = stars["ra"].astype(np.float64) 

720 fgcm_obj["dec"] = stars["dec"].astype(np.float64) 

721 fgcm_obj["obs_arr_index"] = stars["source_cat_index"].astype(np.int64) 

722 fgcm_obj["n_obs"] = stars["nsource"].astype(np.int32) 

723 fgcm_obj["obj_flag"] = np.zeros(len(stars), dtype=np.int32) 

724 

725 # Offset indexes to account for tract merging 

726 fgcm_obj["obs_arr_index"][:] += merge_source_counter 

727 

728 fgcm_objs.append(fgcm_obj) 

729 star_obs_cats.append(star_obs) 

730 

731 merge_source_counter += len(star_obs) 

732 

733 fgcm_obj = vstack(fgcm_objs) 

734 

735 # Set the fgcm_id to a unique 64-bit integer for easier sorting. 

736 fgcm_obj["fgcm_id"][:] = np.arange(len(fgcm_obj)) + 1 

737 

738 if return_observations: 

739 return fgcm_obj, vstack(star_obs_cats) 

740 else: 

741 return fgcm_obj 

742 

743 def _density_downsample(self, fgcm_obj): 

744 """Downsample stars according to density. 

745 

746 Parameters 

747 ---------- 

748 fgcm_obj : `astropy.table.Table` 

749 Catalog of per-star ids and positions. 

750 

751 Returns 

752 ------- 

753 star_ids : `np.ndarray` 

754 Array of downsampled isolated star ids. 

755 """ 

756 if self.config.randomSeed is not None: 

757 rng = np.random.Generator(np.random.MT19937(self.config.randomSeed)) 

758 else: 

759 rng = np.random.Generator(np.random.MT19937()) 

760 

761 ipnest = hpg.angle_to_pixel(self.config.densityCutNside, fgcm_obj["ra"], fgcm_obj["dec"]) 

762 # Use the esutil.stat.histogram function to pull out the histogram of 

763 # grouped pixels, and the rev_indices which describes which inputs 

764 # are grouped together. The fgcm histogram_rev_sorted shim 

765 # ensures that the indices are properly sorted. 

766 hist, rev_indices = histogram_rev_sorted(ipnest) 

767 

768 obj_use = np.ones(len(fgcm_obj), dtype=bool) 

769 

770 (high,) = (hist > self.config.densityCutMaxPerPixel).nonzero() 

771 (ok,) = (hist > 0).nonzero() 

772 self.log.info("There are %d/%d pixels with high stellar density.", high.size, ok.size) 

773 for i in range(high.size): 

774 # The pix_indices are the indices of every star in the pixel. 

775 pix_indices = rev_indices[rev_indices[high[i]]: rev_indices[high[i] + 1]] 

776 # Cut down to the maximum number of stars in the pixel. 

777 cut = rng.choice( 

778 pix_indices, 

779 size=pix_indices.size - self.config.densityCutMaxPerPixel, 

780 replace=False, 

781 ) 

782 obj_use[cut] = False 

783 

784 return np.asarray(fgcm_obj["isolated_star_id"][obj_use]) 

785 

786 def _mark_reserve_stars(self, fgcm_obj): 

787 """Run the star reservation task to flag reserved stars. 

788 

789 Parameters 

790 ---------- 

791 fgcm_obj : `astropy.table.Table` 

792 Catalog of per-star ids and positions. 

793 """ 

794 reserved = self.reserve_selection.run( 

795 len(fgcm_obj), 

796 extra=','.join(self.config.requiredBands), 

797 ) 

798 fgcm_obj["obj_flag"][reserved] |= objFlagDict["RESERVED"] 

799 

800 def _associate_reference_stars(self, lookup_table_handle, visit_cat, fgcm_obj): 

801 """Associate fgcm isolated stars with reference stars. 

802 

803 Parameters 

804 ---------- 

805 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`, optional 

806 Data reference to fgcm look-up table (used if matching reference stars). 

807 visit_cat : `lsst.afw.table.ExposureCatalog` 

808 Catalog of per-visit data. 

809 fgcm_obj : `astropy.table.Table` 

810 Catalog of per-star ids and positions 

811 

812 Returns 

813 ------- 

814 ref_cat : `astropy.table.Table` 

815 Catalog of reference stars matched to fgcm stars. 

816 """ 

817 # Figure out the correct bands/filters that we need based on the data. 

818 lut_cat = lookup_table_handle.get() 

819 

820 std_filter_dict = {filter_name: std_filter for (filter_name, std_filter) in 

821 zip(lut_cat[0]["physicalFilters"].split(","), 

822 lut_cat[0]["stdPhysicalFilters"].split(","))} 

823 std_lambda_dict = {std_filter: std_lambda for (std_filter, std_lambda) in 

824 zip(lut_cat[0]["stdPhysicalFilters"].split(","), 

825 lut_cat[0]["lambdaStdFilter"])} 

826 del lut_cat 

827 

828 reference_filter_names = self._getReferenceFilterNames( 

829 visit_cat, 

830 std_filter_dict, 

831 std_lambda_dict, 

832 ) 

833 self.log.info("Using the following reference filters: %s", (", ".join(reference_filter_names))) 

834 

835 # Split into healpix pixels for more efficient matching.