Coverage for python / lsst / meas / algorithms / loadReferenceObjects.py: 10%

1# This file is part of meas_algorithms. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

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

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

7# for details of code ownership. 

8# 

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

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

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

12# (at your option) any later version. 

13# 

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

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

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

17# GNU General Public License for more details. 

18# 

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

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

21 

22__all__ = ["getRefFluxField", "getRefFluxKeys", "LoadReferenceObjectsConfig", 

23 "ReferenceObjectLoader"] 

24 

25import logging 

26 

27import astropy.time 

28import astropy.units 

29import numpy 

30 

31import lsst.geom as geom 

32import lsst.afw.table as afwTable 

33import lsst.pex.config as pexConfig 

34import lsst.pipe.base as pipeBase 

35from lsst import sphgeom 

36from lsst.daf.base import PropertyList 

37 

38from .convertReferenceCatalog import LATEST_FORMAT_VERSION 

39 

40 

41def getFormatVersionFromRefCat(refCat): 

42 """"Return the format version stored in a reference catalog header. 

43 

44 Parameters 

45 ---------- 

46 refCat : `lsst.afw.table.SimpleCatalog` 

47 Reference catalog to inspect. 

48 

49 Returns 

50 ------- 

51 version : `int` 

52 Format version integer. 

53 

54 Raises 

55 ------ 

56 ValueError 

57 Raised if the catalog is version 0, has no metadata, or does not 

58 include a "REFCAT_FORMAT_VERSION" key. 

59 """ 

60 errMsg = "Version 0 refcats are no longer supported: refcat fluxes must have nJy units." 

61 md = refCat.getMetadata() 

62 if md is None: 

63 raise ValueError(f"No metadata found in refcat header. {errMsg}") 

64 

65 try: 

66 version = md.getScalar("REFCAT_FORMAT_VERSION") 

67 if version == 0: 

68 raise ValueError(errMsg) 

69 else: 

70 return version 

71 except KeyError: 

72 raise ValueError(f"No version number found in refcat header metadata. {errMsg}") 

73 

74 

75class _FilterCatalog: 

76 """This is a private helper class which filters catalogs by 

77 row based on the row being inside the region used to initialize 

78 the class. 

79 

80 Parameters 

81 ---------- 

82 region : `lsst.sphgeom.Region` 

83 The spatial region which all objects should lie within 

84 """ 

85 def __init__(self, region): 

86 self.region = region 

87 

88 def __call__(self, refCat, catRegion): 

89 """This call method on an instance of this class takes in a reference 

90 catalog, and the region from which the catalog was generated. 

91 

92 If the catalog region is entirely contained within the region used to 

93 initialize this class, then all the entries in the catalog must be 

94 within the region and so the whole catalog is returned. 

95 

96 If the catalog region is not entirely contained, then the location for 

97 each record is tested against the region used to initialize the class. 

98 Records which fall inside this region are added to a new catalog, and 

99 this catalog is then returned. 

100 

101 Parameters 

102 --------- 

103 refCat : `lsst.afw.table.SourceCatalog` 

104 SourceCatalog to be filtered. 

105 catRegion : `lsst.sphgeom.Region` 

106 Region in which the catalog was created 

107 """ 

108 if catRegion.isWithin(self.region): 

109 # no filtering needed, region completely contains refcat 

110 return refCat 

111 

112 coordKey = refCat.getCoordKey() 

113 inside = self.region.contains(lon=refCat[coordKey.getRa()], lat=refCat[coordKey.getDec()]) 

114 filteredRefCat = refCat[inside] 

115 

116 return filteredRefCat 

117 

118 

119class LoadReferenceObjectsConfig(pexConfig.Config): 

120 pixelMargin = pexConfig.RangeField( 

121 doc="Padding to add to 4 all edges of the bounding box (pixels)", 

122 dtype=int, 

123 default=250, 

124 min=0, 

125 ) 

126 anyFilterMapsToThis = pexConfig.Field( 

127 doc=("Always use this reference catalog filter, no matter whether or what filter name is " 

128 "supplied to the loader. Effectively a trivial filterMap: map all filter names to this filter." 

129 " This can be set for purely-astrometric catalogs (e.g. Gaia DR2) where there is only one " 

130 "reasonable choice for every camera filter->refcat mapping, but not for refcats used for " 

131 "photometry, which need a filterMap and/or colorterms/transmission corrections."), 

132 dtype=str, 

133 default=None, 

134 optional=True 

135 ) 

136 filterMap = pexConfig.DictField( 

137 doc=("Mapping of camera filter name: reference catalog filter name; " 

138 "each reference filter must exist in the refcat." 

139 " Note that this does not perform any bandpass corrections: it is just a lookup."), 

140 keytype=str, 

141 itemtype=str, 

142 default={}, 

143 ) 

144 requireProperMotion = pexConfig.Field( 

145 doc="Require that the fields needed to correct proper motion " 

146 "(epoch, pm_ra and pm_dec) are present?", 

147 dtype=bool, 

148 default=False, 

149 ) 

150 maxRefObjects = pexConfig.Field( 

151 doc="Maximum number of reference objects to send to the matcher. Setting " 

152 "this to a reasonable value may be desirable for memory reasons " 

153 "(particularly in very crowded field).", 

154 dtype=int, 

155 default=None, 

156 optional=True, 

157 ) 

158 minRefMag = pexConfig.Field( 

159 doc="Minimum (i.e. brightest) magnitude for reference catalog (the brightest " 

160 "sources are typically saturated in the images, so may as well remove " 

161 "them from the reference catalog).", 

162 dtype=float, 

163 default=None, 

164 optional=True, 

165 ) 

166 

167 def validate(self): 

168 super().validate() 

169 if self.filterMap != {} and self.anyFilterMapsToThis is not None: 

170 msg = "`filterMap` and `anyFilterMapsToThis` are mutually exclusive" 

171 raise pexConfig.FieldValidationError(LoadReferenceObjectsConfig.anyFilterMapsToThis, 

172 self, msg) 

173 

174 

175class ReferenceObjectLoader: 

176 """This class facilitates loading reference catalogs. 

177 

178 The QuantumGraph generation will create a list of datasets that may 

179 possibly overlap a given region. These datasets are then used to construct 

180 an instance of this class. The class instance should then be passed into 

181 a task which needs reference catalogs. These tasks should then determine 

182 the exact region of the sky reference catalogs will be loaded for, and 

183 call a corresponding method to load the reference objects. 

184 

185 Parameters 

186 ---------- 

187 dataIds : iterable of `lsst.daf.butler.DataCoordinate` 

188 An iterable object of data IDs that point to reference catalogs. 

189 refCats : iterable of `lsst.daf.butler.DeferredDatasetHandle` 

190 Handles to load refCats on demand. 

191 name : `str`, optional 

192 The name of the refcat that this object will load. This name is used 

193 for applying colorterms, for example. 

194 config : `LoadReferenceObjectsConfig` 

195 Configuration of this reference loader. 

196 log : `lsst.log.Log`, `logging.Logger` or `None`, optional 

197 Logger object used to write out messages. If `None` a default 

198 logger will be used. 

199 """ 

200 ConfigClass = LoadReferenceObjectsConfig 

201 

202 def __init__(self, dataIds, refCats, name=None, log=None, config=None): 

203 if config is None: 

204 config = self.ConfigClass() 

205 self.config = config 

206 self.dataIds = dataIds 

207 self.refCats = list(refCats) 

208 self.name = name 

209 self.log = log or logging.getLogger(__name__).getChild("ReferenceObjectLoader") 

210 

211 def applyProperMotions(self, catalog, epoch): 

212 """Apply proper motion correction to a reference catalog. 

213 

214 Adjust position and position error in the ``catalog`` 

215 for proper motion to the specified ``epoch``, 

216 modifying the catalog in place. 

217 

218 Parameters 

219 ---------- 

220 catalog : `lsst.afw.table.SimpleCatalog` 

221 Catalog of positions, containing at least these fields: 

222 

223 - Coordinates, retrieved by the table's coordinate key. 

224 - ``coord_raErr`` : Error in Right Ascension (rad). 

225 - ``coord_decErr`` : Error in Declination (rad). 

226 - ``pm_ra`` : Proper motion in Right Ascension (rad/yr, 

227 East positive) 

228 - ``pm_raErr`` : Error in ``pm_ra`` (rad/yr), optional. 

229 - ``pm_dec`` : Proper motion in Declination (rad/yr, 

230 North positive) 

231 - ``pm_decErr`` : Error in ``pm_dec`` (rad/yr), optional. 

232 - ``epoch`` : Mean epoch of object (an astropy.time.Time) 

233 epoch : `astropy.time.Time` 

234 Epoch to which to correct proper motion. 

235 If None, do not apply PM corrections or raise if 

236 ``config.requireProperMotion`` is True. 

237 

238 Raises 

239 ------ 

240 RuntimeError 

241 Raised if ``config.requireProperMotion`` is set but we cannot 

242 apply the proper motion correction for some reason. 

243 """ 

244 if epoch is None: 

245 if self.config.requireProperMotion: 

246 raise RuntimeError("requireProperMotion=True but epoch not provided to loader.") 

247 else: 

248 self.log.debug("No epoch provided: not applying proper motion corrections to refcat.") 

249 return 

250 

251 # Warn/raise for a catalog in an incorrect format, if epoch was 

252 # specified. 

253 if "pm_ra" in catalog.schema: 

254 pm_ra_radians = False 

255 field = catalog.schema["pm_ra"].asField() 

256 if field.getTypeString() == "Angle" or field.getUnits() == "rad": 

257 pm_ra_radians = True 

258 

259 if self.config.requireProperMotion and not pm_ra_radians: 

260 raise RuntimeError( 

261 "requireProperMotion=True but refcat pm_ra field is not an Angle or radians.", 

262 ) 

263 elif not pm_ra_radians: 

264 self.log.warning( 

265 "Reference catalog pm_ra field is not an Angle or radians; cannot apply proper motion.", 

266 ) 

267 return 

268 

269 if ("epoch" not in catalog.schema or "pm_ra" not in catalog.schema): 

270 if self.config.requireProperMotion: 

271 raise RuntimeError("requireProperMotion=True but PM data not available from catalog.") 

272 else: 

273 self.log.warning("Proper motion correction not available for this reference catalog.") 

274 return 

275 

276 applyProperMotionsImpl(self.log, catalog, epoch) 

277 

278 @staticmethod 

279 def _remapReferenceCatalogSchema(refCat, *, anyFilterMapsToThis=None, 

280 filterMap=None, centroids=False): 

281 """This function takes in a reference catalog and returns a new catalog 

282 with additional columns defined from the remaining function arguments. 

283 

284 Parameters 

285 ---------- 

286 refCat : `lsst.afw.table.SimpleCatalog` 

287 Reference catalog to map to new catalog 

288 anyFilterMapsToThis : `str`, optional 

289 Always use this reference catalog filter. 

290 Mutually exclusive with `filterMap` 

291 filterMap : `dict` [`str`,`str`], optional 

292 Mapping of camera filter name: reference catalog filter name. 

293 centroids : `bool`, optional 

294 Add centroid fields to the loaded Schema. ``loadPixelBox`` expects 

295 these fields to exist. 

296 

297 Returns 

298 ------- 

299 expandedCat : `lsst.afw.table.SimpleCatalog` 

300 Deep copy of input reference catalog with additional columns added 

301 """ 

302 if anyFilterMapsToThis or filterMap: 

303 ReferenceObjectLoader._addFluxAliases(refCat.schema, anyFilterMapsToThis, filterMap) 

304 

305 mapper = afwTable.SchemaMapper(refCat.schema, True) 

306 mapper.addMinimalSchema(refCat.schema, True) 

307 mapper.editOutputSchema().disconnectAliases() 

308 

309 if centroids: 

310 # Add and initialize centroid and hasCentroid fields (these 

311 # are added after loading to avoid wasting space in the saved 

312 # catalogs). The new fields are automatically initialized to 

313 # (nan, nan) and False so no need to set them explicitly. 

314 mapper.editOutputSchema().addField("centroid_x", type=float, doReplace=True) 

315 mapper.editOutputSchema().addField("centroid_y", type=float, doReplace=True) 

316 mapper.editOutputSchema().addField("hasCentroid", type="Flag", doReplace=True) 

317 mapper.editOutputSchema().getAliasMap().set("slot_Centroid", "centroid") 

318 

319 expandedCat = afwTable.SimpleCatalog(mapper.getOutputSchema()) 

320 expandedCat.setMetadata(refCat.getMetadata()) 

321 expandedCat.extend(refCat, mapper=mapper) 

322 

323 return expandedCat 

324 

325 @staticmethod 

326 def _addFluxAliases(schema, anyFilterMapsToThis=None, filterMap=None): 

327 """Add aliases for camera filter fluxes to the schema. 

328 

329 For each camFilter: refFilter in filterMap, adds these aliases: 

330 <camFilter>_camFlux: <refFilter>_flux 

331 <camFilter>_camFluxErr: <refFilter>_fluxErr, if the latter exists 

332 or sets `anyFilterMapsToThis` in the schema. 

333 

334 Parameters 

335 ---------- 

336 schema : `lsst.afw.table.Schema` 

337 Schema for reference catalog. 

338 anyFilterMapsToThis : `str`, optional 

339 Always use this reference catalog filter. 

340 Mutually exclusive with `filterMap`. 

341 filterMap : `dict` [`str`,`str`], optional 

342 Mapping of camera filter name: reference catalog filter name. 

343 Mutually exclusive with `anyFilterMapsToThis`. 

344 

345 Raises 

346 ------ 

347 RuntimeError 

348 Raised if any required reference flux field is missing from the 

349 schema. 

350 """ 

351 # Fail on any truthy value for either of these. 

352 if anyFilterMapsToThis and filterMap: 

353 raise ValueError("anyFilterMapsToThis and filterMap are mutually exclusive!") 

354 

355 aliasMap = schema.getAliasMap() 

356 

357 if anyFilterMapsToThis is not None: 

358 refFluxName = anyFilterMapsToThis + "_flux" 

359 if refFluxName not in schema: 

360 msg = f"Unknown reference filter for anyFilterMapsToThis='{refFluxName}'" 

361 raise RuntimeError(msg) 

362 aliasMap.set("anyFilterMapsToThis", refFluxName) 

363 return # this is mutually exclusive with filterMap 

364 

365 def addAliasesForOneFilter(filterName, refFilterName): 

366 """Add aliases for a single filter 

367 

368 Parameters 

369 ---------- 

370 filterName : `str` (optional) 

371 Camera filter name. The resulting alias name is 

372 <filterName>_camFlux 

373 refFilterName : `str` 

374 Reference catalog filter name; the field 

375 <refFilterName>_flux must exist. 

376 """ 

377 camFluxName = filterName + "_camFlux" 

378 refFluxName = refFilterName + "_flux" 

379 if refFluxName not in schema: 

380 raise RuntimeError("Unknown reference filter %s" % (refFluxName,)) 

381 aliasMap.set(camFluxName, refFluxName) 

382 refFluxErrName = refFluxName + "Err" 

383 if refFluxErrName in schema: 

384 camFluxErrName = camFluxName + "Err" 

385 aliasMap.set(camFluxErrName, refFluxErrName) 

386 

387 if filterMap is not None: 

388 for filterName, refFilterName in filterMap.items(): 

389 addAliasesForOneFilter(filterName, refFilterName) 

390 

391 @staticmethod 

392 def _makeBoxRegion(BBox, wcs, BBoxPadding): 

393 outerLocalBBox = geom.Box2D(BBox) 

394 innerLocalBBox = geom.Box2D(BBox) 

395 

396 # Grow the bounding box to allow for effects not fully captured by the 

397 # wcs provided (which represents the current best-guess wcs solution 

398 # associated with the dataset for which the calibration is to be 

399 # computed using the loaded and trimmed reference catalog being defined 

400 # here). These effects could include pointing errors and/or an 

401 # insufficient optical distorition model for the instrument. The idea 

402 # is to ensure the spherical geometric region created contains the 

403 # entire region covered by the bbox. 

404 # Also create an inner region that is sure to be inside the bbox. 

405 outerLocalBBox.grow(BBoxPadding) 

406 innerLocalBBox.grow(-1*BBoxPadding) 

407 

408 # Handle the case where the inner bounding box shrank to a zero sized 

409 # region (which will be the case if the shrunken size of either 

410 # dimension is less than or equal to zero). In this case, the inner 

411 # bounding box is set to the original input bounding box. This is 

412 # probably not the best way to handle an empty inner bounding box, but 

413 # it is what the calling code currently expects. 

414 if innerLocalBBox.getDimensions() == geom.Extent2D(0, 0): 

415 innerLocalBBox = geom.Box2D(BBox) 

416 

417 # Convert the corners of the bounding boxes to sky coordinates. 

418 innerBoxCorners = innerLocalBBox.getCorners() 

419 innerSphCorners = [wcs.pixelToSky(corner).getVector() for corner in innerBoxCorners] 

420 innerSkyRegion = sphgeom.ConvexPolygon(innerSphCorners) 

421 

422 outerBoxCorners = outerLocalBBox.getCorners() 

423 outerSphCorners = [wcs.pixelToSky(corner).getVector() for corner in outerBoxCorners] 

424 outerSkyRegion = sphgeom.ConvexPolygon(outerSphCorners) 

425 

426 return innerSkyRegion, outerSkyRegion, innerSphCorners, outerSphCorners 

427 

428 @staticmethod 

429 def _calculateCircle(bbox, wcs, pixelMargin): 

430 """Compute on-sky center and radius of search region. 

431 

432 Parameters 

433 ---------- 

434 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D` 

435 Pixel bounding box. 

436 wcs : `lsst.afw.geom.SkyWcs` 

437 WCS; used to convert pixel positions to sky coordinates. 

438 pixelMargin : `int` 

439 Padding to add to 4 all edges of the bounding box (pixels). 

440 

441 Returns 

442 ------- 

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

444 A Struct containing: 

445 

446 - coord : `lsst.geom.SpherePoint` 

447 ICRS center of the search region. 

448 - radius : `lsst.geom.Angle` 

449 Radius of the search region. 

450 - bbox : `lsst.geom.Box2D` 

451 Bounding box used to compute the circle. 

452 """ 

453 bbox = geom.Box2D(bbox) # we modify the box, so use a copy 

454 bbox.grow(pixelMargin) 

455 coord = wcs.pixelToSky(bbox.getCenter()) 

456 radius = max(coord.separation(wcs.pixelToSky(pp)) for pp in bbox.getCorners()) 

457 return pipeBase.Struct(coord=coord, radius=radius, bbox=bbox) 

458 

459 @staticmethod 

460 def getMetadataCircle(coord, radius, filterName, epoch=None): 

461 """Return metadata about the loaded reference catalog, in an on-sky 

462 circle. 

463 

464 This metadata is used for reloading the catalog (e.g. for 

465 reconstituting a normalized match list). 

466 

467 Parameters 

468 ---------- 

469 coord : `lsst.geom.SpherePoint` 

470 ICRS center of the search region. 

471 radius : `lsst.geom.Angle` 

472 Radius of the search region. 

473 filterName : `str` 

474 Name of the camera filter. 

475 epoch : `astropy.time.Time` or `None`, optional 

476 Epoch that proper motion and parallax were corrected to, or `None` 

477 if no such corrections were applied. 

478 

479 Returns 

480 ------- 

481 md : `lsst.daf.base.PropertyList` 

482 Metadata about the catalog. 

483 """ 

484 md = PropertyList() 

485 md.add('RA', coord.getRa().asDegrees(), 'field center in degrees') 

486 md.add('DEC', coord.getDec().asDegrees(), 'field center in degrees') 

487 md.add('RADIUS', radius.asDegrees(), 'field radius in degrees, minimum') 

488 # Version 1: Initial version 

489 # Version 2: JEPOCH for TAI Julian Epoch year of PM/parallax correction 

490 md.add('SMATCHV', 2, 'SourceMatchVector version number') 

491 md.add('FILTER', filterName, 'camera filter name for photometric data') 

492 md.add('TIMESYS', "TAI", "time scale of time keywords") 

493 md.add('JEPOCH', None if epoch is None else epoch.tai.jyear, 

494 'Julian epoch (TAI Julian Epoch year) for catalog') 

495 return md 

496 

497 def getMetadataBox(self, bbox, wcs, filterName, epoch=None, 

498 bboxToSpherePadding=100): 

499 """Return metadata about the loaded reference catalog, in an 

500 on-detector box. 

501 

502 This metadata is used for reloading the catalog (e.g., for 

503 reconstituting a normalised match list). 

504 

505 Parameters 

506 ---------- 

507 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D` 

508 Bounding box for the pixels. 

509 wcs : `lsst.afw.geom.SkyWcs` 

510 The WCS object associated with ``bbox``. 

511 filterName : `str` 

512 Name of the camera filter. 

513 epoch : `astropy.time.Time` or `None`, optional 

514 Epoch that proper motion and parallax were corrected to, or `None` 

515 if no such corrections were applied. 

516 bboxToSpherePadding : `int`, optional 

517 Padding in pixels to account for translating a set of corners into 

518 a spherical (convex) boundary that is certain to encompass the 

519 enitre area covered by the box. 

520 

521 Returns 

522 ------- 

523 md : `lsst.daf.base.PropertyList` 

524 The metadata detailing the search parameters used for this 

525 dataset. 

526 """ 

527 circle = self._calculateCircle(bbox, wcs, self.config.pixelMargin) 

528 md = self.getMetadataCircle(circle.coord, circle.radius, filterName, epoch=epoch) 

529 

530 paddedBbox = circle.bbox 

531 _, _, innerCorners, outerCorners = self._makeBoxRegion(paddedBbox, wcs, bboxToSpherePadding) 

532 for box, corners in zip(("INNER", "OUTER"), (innerCorners, outerCorners)): 

533 for (name, corner) in zip(("UPPER_LEFT", "UPPER_RIGHT", "LOWER_LEFT", "LOWER_RIGHT"), 

534 corners): 

535 md.add(f"{box}_{name}_RA", geom.SpherePoint(corner).getRa().asDegrees(), f"{box}_corner") 

536 md.add(f"{box}_{name}_DEC", geom.SpherePoint(corner).getDec().asDegrees(), f"{box}_corner") 

537 return md 

538 

539 def loadPixelBox(self, bbox, wcs, filterName, epoch=None, 

540 bboxToSpherePadding=100): 

541 """Load reference objects that are within a pixel-based rectangular 

542 region. 

543 

544 This algorithm works by creating a spherical box whose corners 

545 correspond to the WCS converted corners of the input bounding box 

546 (possibly padded). It then defines a filtering function which looks at 

547 the pixel position of the reference objects and accepts only those that 

548 lie within the specified bounding box. 

549 

550 The spherical box region and filtering function are passed to the 

551 generic loadRegion method which loads and filters the reference objects 

552 from the datastore and returns a single catalog containing the filtered 

553 set of reference objects. 

554 

555 Parameters 

556 ---------- 

557 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D` 

558 Box which bounds a region in pixel space. 

559 wcs : `lsst.afw.geom.SkyWcs` 

560 Wcs object defining the pixel to sky (and inverse) transform for 

561 the supplied ``bbox``. 

562 filterName : `str` 

563 Name of camera filter. 

564 epoch : `astropy.time.Time` or `None`, optional 

565 Epoch to which to correct proper motion and parallax, or `None` 

566 to not apply such corrections. 

567 bboxToSpherePadding : `int`, optional 

568 Padding to account for translating a set of corners into a 

569 spherical (convex) boundary that is certain to encompase the 

570 enitre area covered by the box. 

571 

572 Returns 

573 ------- 

574 output : `lsst.pipe.base.Struct` 

575 Results struct with attributes: 

576 

577 ``refCat`` 

578 Catalog containing reference objects inside the specified 

579 bounding box (padded by self.config.pixelMargin). 

580 ``fluxField`` 

581 Name of the field containing the flux associated with 

582 ``filterName``. 

583 

584 Raises 

585 ------ 

586 RuntimeError 

587 Raised if no reference catalogs could be found for the specified 

588 region. 

589 TypeError 

590 Raised if the loaded reference catalogs do not have matching 

591 schemas. 

592 """ 

593 paddedBbox = geom.Box2D(bbox) 

594 paddedBbox.grow(self.config.pixelMargin) 

595 innerSkyRegion, outerSkyRegion, _, _ = self._makeBoxRegion(paddedBbox, wcs, bboxToSpherePadding) 

596 

597 def _filterFunction(refCat, region): 

598 # Perform an initial "pre filter" step based on the refCat coords 

599 # and the outerSkyRegion created from the self.config.pixelMargin- 

600 # paddedBbox plus an "extra" padding of bboxToSpherePadding and the 

601 # raw wcs. This should ensure a large enough projected area on the 

602 # sky that accounts for any projection/distortion issues, but small 

603 # enough to filter out loaded reference objects that lie well 

604 # beyond the projected detector of interest. This step is required 

605 # due to the very local nature of the wcs available for the 

606 # sky <--> pixel conversions. 

607 preFiltFunc = _FilterCatalog(outerSkyRegion) 

608 refCat = preFiltFunc(refCat, region) 

609 

610 # Add columns to the pre-filtered reference catalog relating their 

611 # coordinates to equivalent pixel positions for the wcs provided 

612 # and use to populate those columns. 

613 refCat = self._remapReferenceCatalogSchema(refCat, centroids=True) 

614 afwTable.updateRefCentroids(wcs, refCat) 

615 # No need to filter the catalog if it is entirely contained in the 

616 # region defined by the inner sky region. 

617 if innerSkyRegion.contains(region): 

618 return refCat 

619 

620 inside = paddedBbox.contains(x=refCat["slot_Centroid_x"], y=refCat["slot_Centroid_y"]) 

621 filteredRefCat = refCat[inside] 

622 

623 return filteredRefCat 

624 return self.loadRegion(outerSkyRegion, filterName, filtFunc=_filterFunction, epoch=epoch) 

625 

626 def loadRegion(self, region, filterName, filtFunc=None, epoch=None, wcsForCentroids=None): 

627 """Load reference objects within a specified region. 

628 

629 This function loads the DataIds used to construct an instance of this 

630 class which intersect or are contained within the specified region. The 

631 reference catalogs which intersect but are not fully contained within 

632 the input region are further filtered by the specified filter function. 

633 This function returns a single source catalog containing all reference 

634 objects inside the specified region. 

635 

636 Parameters 

637 ---------- 

638 region : `lsst.sphgeom.Region` 

639 This can be any type that is derived from `lsst.sphgeom.Region` and 

640 should define the spatial region for which reference objects are to 

641 be loaded. 

642 filtFunc : callable or `None`, optional 

643 This optional parameter should be a callable object that takes a 

644 reference catalog and its corresponding region as parameters, 

645 filters the catalog by some criteria and returns the filtered 

646 reference catalog. If `None`, an internal filter function is used 

647 which filters according to if a reference object falls within the 

648 input region. 

649 filterName : `str` 

650 Name of camera filter. 

651 epoch : `astropy.time.Time` or `None`, optional 

652 Epoch to which to correct proper motion and parallax, or `None` to 

653 not apply such corrections. 

654 wcsForCentroids : `lsst.afw.geom.SkyWcs`, optional 

655 If provided, this WCS will be used to convert the reference catalog 

656 sky coordinates to pixel positions. 

657 

658 Returns 

659 ------- 

660 output : `lsst.pipe.base.Struct` 

661 Results struct with attributes: 

662 

663 ``refCat`` 

664 Catalog containing reference objects which intersect the 

665 input region, filtered by the specified filter function. 

666 ``fluxField`` 

667 Name of the field containing the flux associated with 

668 ``filterName``. 

669 

670 Raises 

671 ------ 

672 RuntimeError 

673 Raised if no reference catalogs could be found for the specified 

674 region. 

675 TypeError 

676 Raised if the loaded reference catalogs do not have matching 

677 schemas. 

678 """ 

679 regionLat = region.getBoundingBox().getLat() 

680 regionLon = region.getBoundingBox().getLon() 

681 self.log.info("Loading reference objects from %s in region bounded by " 

682 "[%.8f, %.8f], [%.8f, %.8f] RA Dec", 

683 self.name, 

684 regionLon.getA().asDegrees(), regionLon.getB().asDegrees(), 

685 regionLat.getA().asDegrees(), regionLat.getB().asDegrees()) 

686 if filtFunc is None: 

687 filtFunc = _FilterCatalog(region) 

688 # Filter out all the regions supplied by the constructor that do not 

689 # overlap. 

690 overlapList = [] 

691 for dataId, refCat in zip(self.dataIds, self.refCats): 

692 # SphGeom supports some objects intersecting others, but is not 

693 # symmetric, try the intersect operation in both directions. 

694 try: 

695 intersects = dataId.region.intersects(region) 

696 except TypeError: 

697 intersects = region.intersects(dataId.region) 

698 

699 if intersects: 

700 overlapList.append((dataId, refCat)) 

701 

702 nOverlap = len(overlapList) 

703 if nOverlap == 0: 

704 raise RuntimeError("No reference tables could be found for input region") 

705 

706 if self.config.maxRefObjects is not None: 

707 maxRefObjectsPerInput = int(self.config.maxRefObjects/nOverlap) 

708 else: 

709 maxRefObjectsPerInput = None 

710 

711 if self.config.anyFilterMapsToThis is not None: 

712 refFluxField = self.config.anyFilterMapsToThis + "_flux" 

713 else: 

714 refFluxField = None 

715 

716 firstCat = overlapList[0][1].get() 

717 # Filter catalog if limits were imposed by the maxRefObjects and/or 

718 # minRefMag config settings. 

719 if (refFluxField is not None 

720 and (maxRefObjectsPerInput is not None or self.config.minRefMag is not None)): 

721 firstCat = filterRefCat(firstCat, refFluxField, maxRefObjectsPerInput, 

722 minRefMag=self.config.minRefMag, log=self.log) 

723 refCat = filtFunc(firstCat, overlapList[0][0].region) 

724 trimmedAmount = len(firstCat) - len(refCat) 

725 

726 # Load in the remaining catalogs 

727 for dataId, inputRefCat in overlapList[1:]: 

728 tmpCat = inputRefCat.get() 

729 

730 if tmpCat.schema != firstCat.schema: 

731 raise TypeError("Reference catalogs have mismatching schemas") 

732 

733 if maxRefObjectsPerInput is not None or self.config.minRefMag is not None: 

734 filteredCat = filterRefCat(tmpCat, refFluxField, maxRefObjectsPerInput, 

735 minRefMag=self.config.minRefMag, log=self.log) 

736 else: 

737 filteredCat = tmpCat 

738 filteredCat = filtFunc(filteredCat, dataId.region) 

739 

740 refCat.extend(filteredCat) 

741 trimmedAmount += len(tmpCat) - len(filteredCat) 

742 

743 self.log.debug("Trimmed %d refCat objects lying outside padded region, leaving %d", 

744 trimmedAmount, len(refCat)) 

745 self.log.info("Loaded %d reference objects", len(refCat)) 

746 

747 # Ensure that the loaded reference catalog is continuous in memory 

748 if not refCat.isContiguous(): 

749 refCat = refCat.copy(deep=True) 

750 

751 self.applyProperMotions(refCat, epoch) 

752 

753 if wcsForCentroids is not None: 

754 expandedCat = self._remapReferenceCatalogSchema( 

755 refCat, 

756 anyFilterMapsToThis=self.config.anyFilterMapsToThis, 

757 filterMap=self.config.filterMap, 

758 centroids=True, 

759 ) 

760 afwTable.updateRefCentroids(wcsForCentroids, expandedCat) 

761 else: 

762 expandedCat = self._remapReferenceCatalogSchema( 

763 refCat, 

764 anyFilterMapsToThis=self.config.anyFilterMapsToThis, 

765 filterMap=self.config.filterMap, 

766 ) 

767 

768 # Ensure that the returned reference catalog is continuous in memory 

769 if not expandedCat.isContiguous(): 

770 expandedCat = expandedCat.copy(deep=True) 

771 

772 fluxField = getRefFluxField(expandedCat.schema, filterName) 

773 

774 if expandedCat.schema[fluxField].asField().getUnits() != "nJy": 

775 # if the flux field is not in nJy, check the refcat format version 

776 version = getFormatVersionFromRefCat(refCat) 

777 if version > LATEST_FORMAT_VERSION: 

778 raise ValueError(f"Unsupported refcat format version: {version} > {LATEST_FORMAT_VERSION}.") 

779 

780 return pipeBase.Struct(refCat=expandedCat, fluxField=fluxField) 

781 

782 def loadSkyCircle(self, ctrCoord, radius, filterName, epoch=None): 

783 """Load reference objects that lie within a circular region on the sky. 

784 

785 This method constructs a circular region from an input center and 

786 angular radius, loads reference catalogs which are contained in or 

787 intersect the circle, and filters reference catalogs which intersect 

788 down to objects which lie within the defined circle. 

789 

790 Parameters 

791 ---------- 

792 ctrCoord : `lsst.geom.SpherePoint` 

793 Point defining the center of the circular region. 

794 radius : `lsst.geom.Angle` 

795 Defines the angular radius of the circular region. 

796 filterName : `str` 

797 Name of camera filter. 

798 epoch : `astropy.time.Time` or `None`, optional 

799 Epoch to which to correct proper motion and parallax, or `None` to 

800 not apply such corrections. 

801 

802 Returns 

803 ------- 

804 output : `lsst.pipe.base.Struct` 

805 Results struct with attributes: 

806 

807 ``refCat`` 

808 Catalog containing reference objects inside the specified 

809 search circle. 

810 ``fluxField`` 

811 Name of the field containing the flux associated with 

812 ``filterName``. 

813 """ 

814 centerVector = ctrCoord.getVector() 

815 sphRadius = sphgeom.Angle(radius.asRadians()) 

816 circularRegion = sphgeom.Circle(centerVector, sphRadius) 

817 return self.loadRegion(circularRegion, filterName, epoch=epoch) 

818 

819 def loadSchema(self, filterName): 

820 """Load the schema for the reference catalog. 

821 

822 Parameters 

823 ---------- 

824 filterName : `str` 

825 Name of camera filter. 

826 

827 Returns 

828 ------- 

829 output : `lsst.pipe.base.Struct` 

830 Results struct with attributes: 

831 

832 ``schema`` 

833 Schema of the reference catalogs returned by other 'load' 

834 methods. 

835 ``fluxField`` 

836 Name of the field containing the flux associated with 

837 ``filterName``. 

838 """ 

839 if not self.refCats: 

840 raise RuntimeError("No reference tables could be found.") 

841 # All refcats should have the same schema, so just get the first one. 

842 cat = self.refCats[0].get() 

843 # Replace the original handle with an in-memory one that caches what 

844 # we've already read, since there's a good chance we'll want to read it 

845 # later. 

846 self.refCats[0] = pipeBase.InMemoryDatasetHandle(cat, dataId=self.refCats[0].dataId, copy=False) 

847 emptyCat = type(cat)(cat.table.clone()) 

848 expandedEmptyCat = self._remapReferenceCatalogSchema( 

849 emptyCat, 

850 anyFilterMapsToThis=self.config.anyFilterMapsToThis, 

851 filterMap=self.config.filterMap, 

852 ) 

853 fluxField = getRefFluxField(expandedEmptyCat.schema, filterName) 

854 if expandedEmptyCat.schema[fluxField].asField().getUnits() != "nJy": 

855 # if the flux field is not in nJy, check the refcat format version 

856 version = getFormatVersionFromRefCat(emptyCat) 

857 if version > LATEST_FORMAT_VERSION: 

858 raise ValueError(f"Unsupported refcat format version: {version} > {LATEST_FORMAT_VERSION}.") 

859 return pipeBase.Struct(schema=expandedEmptyCat.schema, fluxField=fluxField) 

860 

861 

862def filterRefCat(refCat, refFluxField, maxRefObjects=None, minRefMag=None, log=None): 

863 """Sub-select a number of reference objects starting from the brightest 

864 and maxing out at the number specified by maxRefObjects. 

865 

866 No further trimming is done if len(refCat) > maxRefObjects after trimming 

867 to minRefMag. 

868 

869 Parameters 

870 ---------- 

871 refCat : `lsst.afw.table.SimpleCatalog` 

872 Catalog of reference objects to trim. 

873 refFluxField : `str` 

874 Field of refCat to use for flux. 

875 maxRefObjects : `int` or `None`, optional 

876 Maximum number of reference objects (i.e. trim refCat down to 

877 this number of objects). 

878 minRefMag : `int` or `None`, optional 

879 Minimum (i.e. brightest) magnitude to include in the reference 

880 catalog. 

881 log : `lsst.log.Log` or `logging.Logger` or `None`, optional 

882 Logger object used to write out messages. If `None`, no messages 

883 will be logged. 

884 

885 Returns 

886 ------- 

887 filteredCat : `lsst.afw.table.SimpleCatalog` 

888 Catalog trimmed to the maximum brightness and/or maximum number set 

889 in the task config from the brightest flux down. 

890 """ 

891 if maxRefObjects is None and minRefMag is None: 

892 if log is not None: 

893 log.debug("No filtering of the reference catalog has been actually requested " 

894 "(i.e maxRefObjects and minRefMag are both `None`) Returning " 

895 "original catalog.") 

896 return refCat 

897 

898 filteredCat = refCat.copy(deep=True) 

899 

900 if minRefMag is not None: 

901 refFlux = filteredCat.get(refFluxField) 

902 refMag = (refFlux*astropy.units.nJy).to_value(astropy.units.ABmag) 

903 if numpy.nanmin(refMag) <= minRefMag: 

904 filteredCat = filteredCat[(refMag > minRefMag)].copy(deep=True) 

905 if log is not None: 

906 log.info("Trimming the loaded reference catalog to %s > %.2f", refFluxField, minRefMag) 

907 

908 if maxRefObjects is not None: 

909 if len(filteredCat) <= maxRefObjects: 

910 if log is not None: 

911 log.debug("Number of reference objects in reference catalog = %d (max is %d). " 

912 "Returning catalog without further filtering.", 

913 len(filteredCat), maxRefObjects) 

914 return filteredCat 

915 if log is not None: 

916 log.info("Trimming number of reference objects in refCat from %d down to %d. ", 

917 len(refCat), maxRefObjects) 

918 if not filteredCat.isContiguous(): 

919 filteredCat = filteredCat.copy(deep=True) 

920 

921 refFlux = filteredCat.get(refFluxField) 

922 sortedRefFlux = refFlux[refFlux.argsort()] 

923 minRefFlux = sortedRefFlux[-(maxRefObjects + 1)] 

924 

925 selected = (filteredCat.get(refFluxField) > minRefFlux) 

926 filteredCat = filteredCat[selected] 

927 

928 if not filteredCat.isContiguous(): 

929 filteredCat = filteredCat.copy(deep=True) 

930 

931 if log is not None: 

932 if len(filteredCat) > 0: 

933 refFlux = filteredCat[refFluxField] 

934 refMag = (refFlux*astropy.units.nJy).to_value(astropy.units.ABmag) 

935 log.info("Reference catalog magnitude range for filter %s after trimming: refMagMin = %.2f; " 

936 "refMagMax = %.2f", refFluxField, numpy.nanmin(refMag), numpy.nanmax(refMag))