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

1# This file is part of meas_algorithms.

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.

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/>.

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

23 "ReferenceObjectLoader"]

25import logging

27import astropy.time

28import astropy.units

29import numpy

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

38from .convertReferenceCatalog import LATEST_FORMAT_VERSION

41def getFormatVersionFromRefCat(refCat):

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

44 Parameters

45 ----------

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

47 Reference catalog to inspect.

49 Returns

50 -------

51 version : `int`

52 Format version integer.

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}")

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}")

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.

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

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.

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.

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

151 def validate(self):

152 super().validate()

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

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

155 raise pexConfig.FieldValidationError(LoadReferenceObjectsConfig.anyFilterMapsToThis,

156 self, msg)

157

158

159class ReferenceObjectLoader:

160 """This class facilitates loading reference catalogs.

161

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

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

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

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

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

167 call a corresponding method to load the reference objects.

168

169 Parameters

170 ----------

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

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

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

174 Handles to load refCats on demand.

175 name : `str`, optional

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

177 for applying colorterms, for example.

178 config : `LoadReferenceObjectsConfig`

179 Configuration of this reference loader.

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

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

182 logger will be used.

183 """

184 ConfigClass = LoadReferenceObjectsConfig

185

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

187 if config is None:

188 config = self.ConfigClass()

189 self.config = config

190 self.dataIds = dataIds

191 self.refCats = list(refCats)

192 self.name = name

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

194

195 def applyProperMotions(self, catalog, epoch):

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

197

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

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

200 modifying the catalog in place.

201

202 Parameters

203 ----------

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

205 Catalog of positions, containing at least these fields:

206

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

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

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

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

211 East positive)

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

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

214 North positive)

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

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

217 epoch : `astropy.time.Time`

218 Epoch to which to correct proper motion.

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

220 ``config.requireProperMotion`` is True.

221

222 Raises

223 ------

224 RuntimeError

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

226 apply the proper motion correction for some reason.

227 """

228 if epoch is None:

229 if self.config.requireProperMotion:

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

231 else:

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

233 return

234

235 # Warn/raise for a catalog in an incorrect format, if epoch was specified.

236 if "pm_ra" in catalog.schema:

237 pm_ra_radians = False

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

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

240 pm_ra_radians = True

241

242 if self.config.requireProperMotion and not pm_ra_radians:

243 raise RuntimeError(

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

245 )

246 elif not pm_ra_radians:

247 self.log.warning(

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

249 )

250 return

251

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

253 if self.config.requireProperMotion:

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

255 else:

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

257 return

258

259 applyProperMotionsImpl(self.log, catalog, epoch)

260

261 @staticmethod

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

263 filterMap=None, centroids=False):

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

265 with additional columns defined from the remaining function arguments.

266

267 Parameters

268 ----------

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

270 Reference catalog to map to new catalog

271 anyFilterMapsToThis : `str`, optional

272 Always use this reference catalog filter.

273 Mutually exclusive with `filterMap`

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

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

276 centroids : `bool`, optional

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

278 these fields to exist.

279

280 Returns

281 -------

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

283 Deep copy of input reference catalog with additional columns added

284 """

285 if anyFilterMapsToThis or filterMap:

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

287

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

289 mapper.addMinimalSchema(refCat.schema, True)

290 mapper.editOutputSchema().disconnectAliases()

291

292 if centroids:

293 # Add and initialize centroid and hasCentroid fields (these are

294 # added after loading to avoid wasting space in the saved catalogs).

295 # The new fields are automatically initialized to (nan, nan) and

296 # False so no need to set them explicitly.

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

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

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

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

301

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

303 expandedCat.setMetadata(refCat.getMetadata())

304 expandedCat.extend(refCat, mapper=mapper)

305

306 return expandedCat

307

308 @staticmethod

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

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

311

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

313 <camFilter>_camFlux: <refFilter>_flux

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

315 or sets `anyFilterMapsToThis` in the schema.

316

317 Parameters

318 ----------

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

320 Schema for reference catalog.

321 anyFilterMapsToThis : `str`, optional

322 Always use this reference catalog filter.

323 Mutually exclusive with `filterMap`.

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

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

326 Mutually exclusive with `anyFilterMapsToThis`.

327

328 Raises

329 ------

330 RuntimeError

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

332 schema.

333 """

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

335 if anyFilterMapsToThis and filterMap:

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

337

338 aliasMap = schema.getAliasMap()

339

340 if anyFilterMapsToThis is not None:

341 refFluxName = anyFilterMapsToThis + "_flux"

342 if refFluxName not in schema:

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

344 raise RuntimeError(msg)

345 aliasMap.set("anyFilterMapsToThis", refFluxName)

346 return # this is mutually exclusive with filterMap

347

348 def addAliasesForOneFilter(filterName, refFilterName):

349 """Add aliases for a single filter

350

351 Parameters

352 ----------

353 filterName : `str` (optional)

354 Camera filter name. The resulting alias name is

355 <filterName>_camFlux

356 refFilterName : `str`

357 Reference catalog filter name; the field

358 <refFilterName>_flux must exist.

359 """

360 camFluxName = filterName + "_camFlux"

361 refFluxName = refFilterName + "_flux"

362 if refFluxName not in schema:

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

364 aliasMap.set(camFluxName, refFluxName)

365 refFluxErrName = refFluxName + "Err"

366 if refFluxErrName in schema:

367 camFluxErrName = camFluxName + "Err"

368 aliasMap.set(camFluxErrName, refFluxErrName)

369

370 if filterMap is not None:

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

372 addAliasesForOneFilter(filterName, refFilterName)

373

374 @staticmethod

375 def _makeBoxRegion(BBox, wcs, BBoxPadding):

376 outerLocalBBox = geom.Box2D(BBox)

377 innerLocalBBox = geom.Box2D(BBox)

378

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

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

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

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

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

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

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

386 # entire region covered by the bbox.

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

388 outerLocalBBox.grow(BBoxPadding)

389 innerLocalBBox.grow(-1*BBoxPadding)

390

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

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

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

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

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

396 # it is what the calling code currently expects.

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

398 innerLocalBBox = geom.Box2D(BBox)

399

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

401 innerBoxCorners = innerLocalBBox.getCorners()

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

403 innerSkyRegion = sphgeom.ConvexPolygon(innerSphCorners)

404

405 outerBoxCorners = outerLocalBBox.getCorners()

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

407 outerSkyRegion = sphgeom.ConvexPolygon(outerSphCorners)

408

409 return innerSkyRegion, outerSkyRegion, innerSphCorners, outerSphCorners

410

411 @staticmethod

412 def _calculateCircle(bbox, wcs, pixelMargin):

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

414

415 Parameters

416 ----------

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

418 Pixel bounding box.

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

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

421 pixelMargin : `int`

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

423

424 Returns

425 -------

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

427 A Struct containing:

428

429 - coord : `lsst.geom.SpherePoint`

430 ICRS center of the search region.

431 - radius : `lsst.geom.Angle`

432 Radius of the search region.

433 - bbox : `lsst.geom.Box2D`

434 Bounding box used to compute the circle.

435 """

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

437 bbox.grow(pixelMargin)

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

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

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

441

442 @staticmethod

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

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

445 circle.

446

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

448 reconstituting a normalized match list).

449

450 Parameters

451 ----------

452 coord : `lsst.geom.SpherePoint`

453 ICRS center of the search region.

454 radius : `lsst.geom.Angle`

455 Radius of the search region.

456 filterName : `str`

457 Name of the camera filter.

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

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

460 if no such corrections were applied.

461

462 Returns

463 -------

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

465 Metadata about the catalog.

466 """

467 md = PropertyList()

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

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

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

471 # Version 1: Initial version

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

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

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

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

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

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

478 return md

479

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

481 bboxToSpherePadding=100):

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

483 on-detector box.

484

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

486 reconstituting a normalised match list).

487

488 Parameters

489 ----------

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

491 Bounding box for the pixels.

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

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

494 filterName : `str`

495 Name of the camera filter.

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

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

498 if no such corrections were applied.

499 bboxToSpherePadding : `int`, optional

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

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

502 enitre area covered by the box.

503

504 Returns

505 -------

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

507 The metadata detailing the search parameters used for this

508 dataset.

509 """

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

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

512

513 paddedBbox = circle.bbox

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

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

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

517 corners):

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

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

520 return md

521

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

523 bboxToSpherePadding=100):

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

525 region.

526

527 This algorithm works by creating a spherical box whose corners

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

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

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

531 lie within the specified bounding box.

532

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

534 generic loadRegion method which loads and filters the reference objects

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

536 set of reference objects.

537

538 Parameters

539 ----------

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

541 Box which bounds a region in pixel space.

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

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

544 the supplied ``bbox``.

545 filterName : `str`

546 Name of camera filter.

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

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

549 to not apply such corrections.

550 bboxToSpherePadding : `int`, optional

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

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

553 enitre area covered by the box.

554

555 Returns

556 -------

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

558 Results struct with attributes:

559

560 ``refCat``

561 Catalog containing reference objects inside the specified

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

563 ``fluxField``

564 Name of the field containing the flux associated with

565 ``filterName``.

566

567 Raises

568 ------

569 RuntimeError

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

571 region.

572 TypeError

573 Raised if the loaded reference catalogs do not have matching

574 schemas.

575 """

576 paddedBbox = geom.Box2D(bbox)

577 paddedBbox.grow(self.config.pixelMargin)

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

579

580 def _filterFunction(refCat, region):

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

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

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

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

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

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

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

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

589 # sky <--> pixel conversions.

590 preFiltFunc = _FilterCatalog(outerSkyRegion)

591 refCat = preFiltFunc(refCat, region)

592

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

594 # coordinates to equivalent pixel positions for the wcs provided

595 # and use to populate those columns.

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

597 afwTable.updateRefCentroids(wcs, refCat)

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

599 # region defined by the inner sky region.

600 if innerSkyRegion.contains(region):

601 return refCat

602

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

604 filteredRefCat = refCat[inside]

605

606 return filteredRefCat

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

608

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

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

611

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

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

614 reference catalogs which intersect but are not fully contained within

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

616 This function returns a single source catalog containing all reference

617 objects inside the specified region.

618

619 Parameters

620 ----------

621 region : `lsst.sphgeom.Region`

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

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

624 be loaded.

625 filtFunc : callable or `None`, optional

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

627 reference catalog and its corresponding region as parameters,

628 filters the catalog by some criteria and returns the filtered

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

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

631 input region.

632 filterName : `str`

633 Name of camera filter.

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

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

636 not apply such corrections.

637

638 Returns

639 -------

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

641 Results struct with attributes:

642

643 ``refCat``

644 Catalog containing reference objects which intersect the

645 input region, filtered by the specified filter function.

646 ``fluxField``

647 Name of the field containing the flux associated with

648 ``filterName``.

649

650 Raises

651 ------

652 RuntimeError

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

654 region.

655 TypeError

656 Raised if the loaded reference catalogs do not have matching

657 schemas.

658 """

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

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

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

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

663 self.name,

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

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

666 if filtFunc is None:

667 filtFunc = _FilterCatalog(region)

668 # filter out all the regions supplied by the constructor that do not overlap

669 overlapList = []

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

671 # SphGeom supports some objects intersecting others, but is not symmetric,

672 # try the intersect operation in both directions

673 try:

674 intersects = dataId.region.intersects(region)

675 except TypeError:

676 intersects = region.intersects(dataId.region)

677

678 if intersects:

679 overlapList.append((dataId, refCat))

680

681 if len(overlapList) == 0:

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

683

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

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

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

687

688 # Load in the remaining catalogs

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

690 tmpCat = inputRefCat.get()

691

692 if tmpCat.schema != firstCat.schema:

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

694

695 filteredCat = filtFunc(tmpCat, dataId.region)

696 refCat.extend(filteredCat)

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

698

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

700 trimmedAmount, len(refCat))

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

702

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

704 if not refCat.isContiguous():

705 refCat = refCat.copy(deep=True)

706

707 self.applyProperMotions(refCat, epoch)

708

709 expandedCat = self._remapReferenceCatalogSchema(refCat,

710 anyFilterMapsToThis=self.config.anyFilterMapsToThis,

711 filterMap=self.config.filterMap)

712

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

714 if not expandedCat.isContiguous():

715 expandedCat = expandedCat.copy(deep=True)

716

717 fluxField = getRefFluxField(expandedCat.schema, filterName)

718

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

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

721 version = getFormatVersionFromRefCat(refCat)

722 if version > LATEST_FORMAT_VERSION:

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

724

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

726

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

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

729

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

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

732 intersect the circle, and filters reference catalogs which intersect

733 down to objects which lie within the defined circle.

734

735 Parameters

736 ----------

737 ctrCoord : `lsst.geom.SpherePoint`

738 Point defining the center of the circular region.

739 radius : `lsst.geom.Angle`

740 Defines the angular radius of the circular region.

741 filterName : `str`

742 Name of camera filter.

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

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

745 not apply such corrections.

746

747 Returns

748 -------

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

750 Results struct with attributes:

751

752 ``refCat``

753 Catalog containing reference objects inside the specified

754 search circle.

755 ``fluxField``

756 Name of the field containing the flux associated with

757 ``filterName``.

758 """

759 centerVector = ctrCoord.getVector()

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

761 circularRegion = sphgeom.Circle(centerVector, sphRadius)

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

763

764 def loadSchema(self, filterName):

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

766

767 Parameters

768 ----------

769 filterName : `str`

770 Name of camera filter.

771

772 Returns

773 -------

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

775 Results struct with attributes:

776

777 ``schema``

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

779 methods.

780 ``fluxField``

781 Name of the field containing the flux associated with

782 ``filterName``.

783 """

784 if not self.refCats:

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

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

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

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

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

790 # later.

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

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

793 expandedEmptyCat = self._remapReferenceCatalogSchema(

794 emptyCat,

795 anyFilterMapsToThis=self.config.anyFilterMapsToThis,

796 filterMap=self.config.filterMap,

797 )

798 fluxField = getRefFluxField(expandedEmptyCat.schema, filterName)

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

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

801 version = getFormatVersionFromRefCat(emptyCat)

802 if version > LATEST_FORMAT_VERSION:

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

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

805

806

807def getRefFluxField(schema, filterName):

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

809

810 Parameters

811 ----------

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

813 Reference catalog schema.

814 filterName : `str`

815 Name of camera filter.

816

817 Returns

818 -------

819 fluxFieldName : `str`

820 Name of flux field.

821

822 Notes

823 -----

824 Return the alias of ``anyFilterMapsToThis``, if present

825 else, return ``*filterName*_camFlux`` if present,

826 else, return ``*filterName*_flux`` if present (camera filter name

827 matches reference filter name), else raise an exception.

828

829 Raises

830 ------

831 RuntimeError

832 Raised if an appropriate field is not found.

833 """

834 if not isinstance(schema, afwTable.Schema):

835 raise RuntimeError("schema=%s is not a schema" % (schema,))

836 try:

837 return schema.getAliasMap().get("anyFilterMapsToThis")

838 except LookupError:

839 pass # try the filterMap next

840

841 fluxFieldList = [filterName + "_camFlux", filterName + "_flux"]

842 for fluxField in fluxFieldList:

843 if fluxField in schema:

844 return fluxField

845

846 raise RuntimeError("Could not find flux field(s) %s" % (", ".join(fluxFieldList)))

847

848

849def getRefFluxKeys(schema, filterName):

850 """Return keys for flux and flux error.

851

852 Parameters

853 ----------

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

855 Reference catalog schema.

856 filterName : `str`

857 Name of camera filter.

858

859 Returns

860 -------

861 keys : `tuple` of (`lsst.afw.table.Key`, `lsst.afw.table.Key`)

862 Two keys:

863

864 - flux key

865 - flux error key, if present, else None

866

867 Raises

868 ------

869 RuntimeError

870 If flux field not found.

871 """

872 fluxField = getRefFluxField(schema, filterName)

873 fluxErrField = fluxField + "Err"

874 fluxKey = schema[fluxField].asKey()

875 try:

876 fluxErrKey = schema[fluxErrField].asKey()

877 except Exception:

878 fluxErrKey = None

879 return (fluxKey, fluxErrKey)

880

881

882def applyProperMotionsImpl(log, catalog, epoch):

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

884

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

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

887 modifying the catalog in place.

888

889 Parameters

890 ----------

891 log : `lsst.log.Log` or `logging.Logger`

892 Log object to write to.

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

894 Catalog of positions, containing:

895

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

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

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

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

900 East positive)

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

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

903 North positive)

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

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

906 epoch : `astropy.time.Time`

907 Epoch to which to correct proper motion.

908 """

909 if "epoch" not in catalog.schema or "pm_ra" not in catalog.schema or "pm_dec" not in catalog.schema:

910 log.warning("Proper motion correction not available from catalog")

911 return

912 if not catalog.isContiguous():

913 raise RuntimeError("Catalog must be contiguous")

914 catEpoch = astropy.time.Time(catalog["epoch"], scale="tai", format="mjd")

915 log.info("Correcting reference catalog for proper motion to %r", epoch)

916 # Use `epoch.tai` to make sure the time difference is in TAI

917 timeDiffsYears = (epoch.tai - catEpoch).to(astropy.units.yr).value

918 coordKey = catalog.table.getCoordKey()

919 # Compute the offset of each object due to proper motion

920 # as components of the arc of a great circle along RA and Dec

921 pmRaRad = catalog["pm_ra"]

922 pmDecRad = catalog["pm_dec"]

923 offsetsRaRad = pmRaRad*timeDiffsYears

924 offsetsDecRad = pmDecRad*timeDiffsYears

925 # Compute the corresponding bearing and arc length of each offset

926 # due to proper motion, and apply the offset.

927 # The factor of 1e6 for computing bearing is intended as

928 # a reasonable scale for typical values of proper motion

929 # in order to avoid large errors for small values of proper motion;

930 # using the offsets is another option, but it can give

931 # needlessly large errors for short duration.

932 offsetBearingsRad = numpy.arctan2(offsetsDecRad*1e6, offsetsRaRad*1e6)

933 offsetAmountsRad = numpy.hypot(offsetsRaRad, offsetsDecRad)

934 for record, bearingRad, amountRad in zip(catalog, offsetBearingsRad, offsetAmountsRad):

935 record.set(coordKey,

936 record.get(coordKey).offset(bearing=bearingRad*geom.radians,

937 amount=amountRad*geom.radians))

938 # TODO DM-36979: this needs to incorporate the full covariance!

939 # Increase error in RA and Dec based on error in proper motion

940 if "coord_raErr" in catalog.schema:

941 catalog["coord_raErr"] = numpy.hypot(catalog["coord_raErr"],

942 catalog["pm_raErr"]*timeDiffsYears)

943 if "coord_decErr" in catalog.schema:

944 catalog["coord_decErr"] = numpy.hypot(catalog["coord_decErr"],

945 catalog["pm_decErr"]*timeDiffsYears)