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): 

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 

655 Returns 

656 ------- 

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

658 Results struct with attributes: 

659 

660 ``refCat`` 

661 Catalog containing reference objects which intersect the 

662 input region, filtered by the specified filter function. 

663 ``fluxField`` 

664 Name of the field containing the flux associated with 

665 ``filterName``. 

666 

667 Raises 

668 ------ 

669 RuntimeError 

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

671 region. 

672 TypeError 

673 Raised if the loaded reference catalogs do not have matching 

674 schemas. 

675 """ 

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

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

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

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

680 self.name, 

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

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

683 if filtFunc is None: 

684 filtFunc = _FilterCatalog(region) 

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

686 # overlap. 

687 overlapList = [] 

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

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

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

691 try: 

692 intersects = dataId.region.intersects(region) 

693 except TypeError: 

694 intersects = region.intersects(dataId.region) 

695 

696 if intersects: 

697 overlapList.append((dataId, refCat)) 

698 

699 nOverlap = len(overlapList) 

700 if nOverlap == 0: 

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

702 

703 if self.config.maxRefObjects is not None: 

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

705 else: 

706 maxRefObjectsPerInput = None 

707 

708 if self.config.anyFilterMapsToThis is not None: 

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

710 else: 

711 refFluxField = None 

712 

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

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

715 # minRefMag config settings. 

716 if (refFluxField is not None 

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

718 firstCat = filterRefCat(firstCat, refFluxField, maxRefObjectsPerInput, 

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

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

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

722 

723 # Load in the remaining catalogs 

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

725 tmpCat = inputRefCat.get() 

726 

727 if tmpCat.schema != firstCat.schema: 

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

729 

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

731 filteredCat = filterRefCat(tmpCat, refFluxField, maxRefObjectsPerInput, 

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

733 else: 

734 filteredCat = tmpCat 

735 filteredCat = filtFunc(filteredCat, dataId.region) 

736 

737 refCat.extend(filteredCat) 

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

739 

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

741 trimmedAmount, len(refCat)) 

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

743 

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

745 if not refCat.isContiguous(): 

746 refCat = refCat.copy(deep=True) 

747 

748 self.applyProperMotions(refCat, epoch) 

749 

750 expandedCat = self._remapReferenceCatalogSchema(refCat, 

751 anyFilterMapsToThis=self.config.anyFilterMapsToThis, 

752 filterMap=self.config.filterMap) 

753 

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

755 if not expandedCat.isContiguous(): 

756 expandedCat = expandedCat.copy(deep=True) 

757 

758 fluxField = getRefFluxField(expandedCat.schema, filterName) 

759 

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

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

762 version = getFormatVersionFromRefCat(refCat) 

763 if version > LATEST_FORMAT_VERSION: 

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

765 

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

767 

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

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

770 

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

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

773 intersect the circle, and filters reference catalogs which intersect 

774 down to objects which lie within the defined circle. 

775 

776 Parameters 

777 ---------- 

778 ctrCoord : `lsst.geom.SpherePoint` 

779 Point defining the center of the circular region. 

780 radius : `lsst.geom.Angle` 

781 Defines the angular radius of the circular region. 

782 filterName : `str` 

783 Name of camera filter. 

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

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

786 not apply such corrections. 

787 

788 Returns 

789 ------- 

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

791 Results struct with attributes: 

792 

793 ``refCat`` 

794 Catalog containing reference objects inside the specified 

795 search circle. 

796 ``fluxField`` 

797 Name of the field containing the flux associated with 

798 ``filterName``. 

799 """ 

800 centerVector = ctrCoord.getVector() 

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

802 circularRegion = sphgeom.Circle(centerVector, sphRadius) 

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

804 

805 def loadSchema(self, filterName): 

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

807 

808 Parameters 

809 ---------- 

810 filterName : `str` 

811 Name of camera filter. 

812 

813 Returns 

814 ------- 

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

816 Results struct with attributes: 

817 

818 ``schema`` 

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

820 methods. 

821 ``fluxField`` 

822 Name of the field containing the flux associated with 

823 ``filterName``. 

824 """ 

825 if not self.refCats: 

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

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

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

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

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

831 # later. 

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

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

834 expandedEmptyCat = self._remapReferenceCatalogSchema( 

835 emptyCat, 

836 anyFilterMapsToThis=self.config.anyFilterMapsToThis, 

837 filterMap=self.config.filterMap, 

838 ) 

839 fluxField = getRefFluxField(expandedEmptyCat.schema, filterName) 

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

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

842 version = getFormatVersionFromRefCat(emptyCat) 

843 if version > LATEST_FORMAT_VERSION: 

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

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

846 

847 

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

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

850 and maxing out at the number specified by maxRefObjects. 

851 

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

853 to minRefMag. 

854 

855 Parameters 

856 ---------- 

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

858 Catalog of reference objects to trim. 

859 refFluxField : `str` 

860 Field of refCat to use for flux. 

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

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

863 this number of objects). 

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

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

866 catalog. 

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

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

869 will be logged. 

870 

871 Returns 

872 ------- 

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

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

875 in the task config from the brightest flux down. 

876 """ 

877 if maxRefObjects is None and minRefMag is None: 

878 if log is not None: 

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

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

881 "original catalog.") 

882 return refCat 

883 

884 filteredCat = refCat.copy(deep=True) 

885 

886 if minRefMag is not None: 

887 refFlux = filteredCat.get(refFluxField) 

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

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

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

891 if log is not None: 

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

893 

894 if maxRefObjects is not None: 

895 if len(filteredCat) <= maxRefObjects: 

896 if log is not None: 

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

898 "Returning catalog without further filtering.", 

899 len(filteredCat), maxRefObjects) 

900 return filteredCat 

901 if log is not None: 

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

903 len(refCat), maxRefObjects) 

904 if not filteredCat.isContiguous(): 

905 filteredCat = filteredCat.copy(deep=True) 

906 

907 refFlux = filteredCat.get(refFluxField) 

908 sortedRefFlux = refFlux[refFlux.argsort()] 

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

910 

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

912 filteredCat = filteredCat[selected] 

913 

914 if not filteredCat.isContiguous(): 

915 filteredCat = filteredCat.copy(deep=True) 

916 

917 if log is not None: 

918 if len(filteredCat) > 0: 

919 refFlux = filteredCat[refFluxField] 

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

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

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

923 else: 

924 log.warning("Length of reference catalog after filtering is 0.") 

925 return filteredCat 

926 

927 

928def getRefFluxField(schema, filterName): 

929 """Get the name of a flux field from a schema. 

930 

931 Parameters 

932 ---------- 

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

934 Reference catalog schema. 

935 filterName : `str` 

936 Name of camera filter. 

937 

938 Returns 

939 ------- 

940 fluxFieldName : `str`