LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
loadReferenceObjects.py
Go to the documentation of this file.
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
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
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 filteredRefCat = type(refCat)(refCat.table)
113 for record in refCat:
114 if self.region.contains(record.getCoord().getVector()):
115 filteredRefCat.append(record)
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
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 = 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 and not isinstance(catalog.schema["pm_ra"].asKey(), afwTable.KeyAngle)):
238 if self.config.requireProperMotion:
239 raise RuntimeError("requireProperMotion=True but refcat pm_ra field is not an Angle.")
240 else:
241 self.log.warning("Reference catalog pm_ra field is not an Angle; cannot apply proper motion.")
242 return
243
244 if ("epoch" not in catalog.schema or "pm_ra" not in catalog.schema):
245 if self.config.requireProperMotion:
246 raise RuntimeError("requireProperMotion=True but PM data not available from catalog.")
247 else:
248 self.log.warning("Proper motion correction not available for this reference catalog.")
249 return
250
251 applyProperMotionsImpl(self.log, catalog, epoch)
252
253 @staticmethod
254 def _remapReferenceCatalogSchema(refCat, *, anyFilterMapsToThis=None,
255 filterMap=None, centroids=False):
256 """This function takes in a reference catalog and returns a new catalog
257 with additional columns defined from the remaining function arguments.
258
259 Parameters
260 ----------
261 refCat : `lsst.afw.table.SimpleCatalog`
262 Reference catalog to map to new catalog
263 anyFilterMapsToThis : `str`, optional
264 Always use this reference catalog filter.
265 Mutually exclusive with `filterMap`
266 filterMap : `dict` [`str`,`str`], optional
267 Mapping of camera filter name: reference catalog filter name.
268 centroids : `bool`, optional
269 Add centroid fields to the loaded Schema. ``loadPixelBox`` expects
270 these fields to exist.
271
272 Returns
273 -------
274 expandedCat : `lsst.afw.table.SimpleCatalog`
275 Deep copy of input reference catalog with additional columns added
276 """
277 if anyFilterMapsToThis or filterMap:
278 ReferenceObjectLoader._addFluxAliases(refCat.schema, anyFilterMapsToThis, filterMap)
279
280 mapper = afwTable.SchemaMapper(refCat.schema, True)
281 mapper.addMinimalSchema(refCat.schema, True)
282 mapper.editOutputSchema().disconnectAliases()
283
284 if centroids:
285 # Add and initialize centroid and hasCentroid fields (these are
286 # added after loading to avoid wasting space in the saved catalogs).
287 # The new fields are automatically initialized to (nan, nan) and
288 # False so no need to set them explicitly.
289 mapper.editOutputSchema().addField("centroid_x", type=float, doReplace=True)
290 mapper.editOutputSchema().addField("centroid_y", type=float, doReplace=True)
291 mapper.editOutputSchema().addField("hasCentroid", type="Flag", doReplace=True)
292 mapper.editOutputSchema().getAliasMap().set("slot_Centroid", "centroid")
293
294 expandedCat = afwTable.SimpleCatalog(mapper.getOutputSchema())
295 expandedCat.setMetadata(refCat.getMetadata())
296 expandedCat.extend(refCat, mapper=mapper)
297
298 return expandedCat
299
300 @staticmethod
301 def _addFluxAliases(schema, anyFilterMapsToThis=None, filterMap=None):
302 """Add aliases for camera filter fluxes to the schema.
303
304 For each camFilter: refFilter in filterMap, adds these aliases:
305 <camFilter>_camFlux: <refFilter>_flux
306 <camFilter>_camFluxErr: <refFilter>_fluxErr, if the latter exists
307 or sets `anyFilterMapsToThis` in the schema.
308
309 Parameters
310 ----------
311 schema : `lsst.afw.table.Schema`
312 Schema for reference catalog.
313 anyFilterMapsToThis : `str`, optional
314 Always use this reference catalog filter.
315 Mutually exclusive with `filterMap`.
316 filterMap : `dict` [`str`,`str`], optional
317 Mapping of camera filter name: reference catalog filter name.
318 Mutually exclusive with `anyFilterMapsToThis`.
319
320 Raises
321 ------
322 RuntimeError
323 Raised if any required reference flux field is missing from the
324 schema.
325 """
326 # Fail on any truthy value for either of these.
327 if anyFilterMapsToThis and filterMap:
328 raise ValueError("anyFilterMapsToThis and filterMap are mutually exclusive!")
329
330 aliasMap = schema.getAliasMap()
331
332 if anyFilterMapsToThis is not None:
333 refFluxName = anyFilterMapsToThis + "_flux"
334 if refFluxName not in schema:
335 msg = f"Unknown reference filter for anyFilterMapsToThis='{refFluxName}'"
336 raise RuntimeError(msg)
337 aliasMap.set("anyFilterMapsToThis", refFluxName)
338 return # this is mutually exclusive with filterMap
339
340 def addAliasesForOneFilter(filterName, refFilterName):
341 """Add aliases for a single filter
342
343 Parameters
344 ----------
345 filterName : `str` (optional)
346 Camera filter name. The resulting alias name is
347 <filterName>_camFlux
348 refFilterName : `str`
349 Reference catalog filter name; the field
350 <refFilterName>_flux must exist.
351 """
352 camFluxName = filterName + "_camFlux"
353 refFluxName = refFilterName + "_flux"
354 if refFluxName not in schema:
355 raise RuntimeError("Unknown reference filter %s" % (refFluxName,))
356 aliasMap.set(camFluxName, refFluxName)
357 refFluxErrName = refFluxName + "Err"
358 if refFluxErrName in schema:
359 camFluxErrName = camFluxName + "Err"
360 aliasMap.set(camFluxErrName, refFluxErrName)
361
362 if filterMap is not None:
363 for filterName, refFilterName in filterMap.items():
364 addAliasesForOneFilter(filterName, refFilterName)
365
366 @staticmethod
367 def _makeBoxRegion(BBox, wcs, BBoxPadding):
368 outerLocalBBox = geom.Box2D(BBox)
369 innerLocalBBox = geom.Box2D(BBox)
370
371 # Grow the bounding box to allow for effects not fully captured by the
372 # wcs provided (which represents the current best-guess wcs solution
373 # associated with the dataset for which the calibration is to be
374 # computed using the loaded and trimmed reference catalog being defined
375 # here). These effects could include pointing errors and/or an
376 # insufficient optical distorition model for the instrument. The idea
377 # is to ensure the spherical geometric region created contains the
378 # entire region covered by the bbox.
379 # Also create an inner region that is sure to be inside the bbox.
380 outerLocalBBox.grow(BBoxPadding)
381 innerLocalBBox.grow(-1*BBoxPadding)
382
383 # Handle the case where the inner bounding box shrank to a zero sized
384 # region (which will be the case if the shrunken size of either
385 # dimension is less than or equal to zero). In this case, the inner
386 # bounding box is set to the original input bounding box. This is
387 # probably not the best way to handle an empty inner bounding box, but
388 # it is what the calling code currently expects.
389 if innerLocalBBox.getDimensions() == geom.Extent2D(0, 0):
390 innerLocalBBox = geom.Box2D(BBox)
391
392 # Convert the corners of the bounding boxes to sky coordinates.
393 innerBoxCorners = innerLocalBBox.getCorners()
394 innerSphCorners = [wcs.pixelToSky(corner).getVector() for corner in innerBoxCorners]
395 innerSkyRegion = sphgeom.ConvexPolygon(innerSphCorners)
396
397 outerBoxCorners = outerLocalBBox.getCorners()
398 outerSphCorners = [wcs.pixelToSky(corner).getVector() for corner in outerBoxCorners]
399 outerSkyRegion = sphgeom.ConvexPolygon(outerSphCorners)
400
401 return innerSkyRegion, outerSkyRegion, innerSphCorners, outerSphCorners
402
403 @staticmethod
404 def _calculateCircle(bbox, wcs, pixelMargin):
405 """Compute on-sky center and radius of search region.
406
407 Parameters
408 ----------
409 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D`
410 Pixel bounding box.
411 wcs : `lsst.afw.geom.SkyWcs`
412 WCS; used to convert pixel positions to sky coordinates.
413 pixelMargin : `int`
414 Padding to add to 4 all edges of the bounding box (pixels).
415
416 Returns
417 -------
418 results : `lsst.pipe.base.Struct`
419 A Struct containing:
420
421 - coord : `lsst.geom.SpherePoint`
422 ICRS center of the search region.
423 - radius : `lsst.geom.Angle`
424 Radius of the search region.
425 - bbox : `lsst.geom.Box2D`
426 Bounding box used to compute the circle.
427 """
428 bbox = geom.Box2D(bbox) # we modify the box, so use a copy
429 bbox.grow(pixelMargin)
430 coord = wcs.pixelToSky(bbox.getCenter())
431 radius = max(coord.separation(wcs.pixelToSky(pp)) for pp in bbox.getCorners())
432 return pipeBase.Struct(coord=coord, radius=radius, bbox=bbox)
433
434 @staticmethod
435 def getMetadataCircle(coord, radius, filterName, epoch=None):
436 """Return metadata about the loaded reference catalog, in an on-sky
437 circle.
438
439 This metadata is used for reloading the catalog (e.g. for
440 reconstituting a normalized match list).
441
442 Parameters
443 ----------
444 coord : `lsst.geom.SpherePoint`
445 ICRS center of the search region.
446 radius : `lsst.geom.Angle`
447 Radius of the search region.
448 filterName : `str`
449 Name of the camera filter.
450 epoch : `astropy.time.Time` or `None`, optional
451 Epoch that proper motion and parallax were corrected to, or `None`
452 if no such corrections were applied.
453
454 Returns
455 -------
456 md : `lsst.daf.base.PropertyList`
457 Metadata about the catalog.
458 """
459 md = PropertyList()
460 md.add('RA', coord.getRa().asDegrees(), 'field center in degrees')
461 md.add('DEC', coord.getDec().asDegrees(), 'field center in degrees')
462 md.add('RADIUS', radius.asDegrees(), 'field radius in degrees, minimum')
463 # Version 1: Initial version
464 # Version 2: JEPOCH for TAI Julian Epoch year of PM/parallax correction
465 md.add('SMATCHV', 2, 'SourceMatchVector version number')
466 md.add('FILTER', filterName, 'camera filter name for photometric data')
467 md.add('TIMESYS', "TAI", "time scale of time keywords")
468 md.add('JEPOCH', None if epoch is None else epoch.tai.jyear,
469 'Julian epoch (TAI Julian Epoch year) for catalog')
470 return md
471
472 def getMetadataBox(self, bbox, wcs, filterName, epoch=None,
473 bboxToSpherePadding=100):
474 """Return metadata about the loaded reference catalog, in an
475 on-detector box.
476
477 This metadata is used for reloading the catalog (e.g., for
478 reconstituting a normalised match list).
479
480 Parameters
481 ----------
482 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D`
483 Bounding box for the pixels.
484 wcs : `lsst.afw.geom.SkyWcs`
485 The WCS object associated with ``bbox``.
486 filterName : `str`
487 Name of the camera filter.
488 epoch : `astropy.time.Time` or `None`, optional
489 Epoch that proper motion and parallax were corrected to, or `None`
490 if no such corrections were applied.
491 bboxToSpherePadding : `int`, optional
492 Padding in pixels to account for translating a set of corners into
493 a spherical (convex) boundary that is certain to encompass the
494 enitre area covered by the box.
495
496 Returns
497 -------
498 md : `lsst.daf.base.PropertyList`
499 The metadata detailing the search parameters used for this
500 dataset.
501 """
502 circle = self._calculateCircle(bbox, wcs, self.config.pixelMargin)
503 md = self.getMetadataCircle(circle.coord, circle.radius, filterName, epoch=epoch)
504
505 paddedBbox = circle.bbox
506 _, _, innerCorners, outerCorners = self._makeBoxRegion(paddedBbox, wcs, bboxToSpherePadding)
507 for box, corners in zip(("INNER", "OUTER"), (innerCorners, outerCorners)):
508 for (name, corner) in zip(("UPPER_LEFT", "UPPER_RIGHT", "LOWER_LEFT", "LOWER_RIGHT"),
509 corners):
510 md.add(f"{box}_{name}_RA", geom.SpherePoint(corner).getRa().asDegrees(), f"{box}_corner")
511 md.add(f"{box}_{name}_DEC", geom.SpherePoint(corner).getDec().asDegrees(), f"{box}_corner")
512 return md
513
514 def loadPixelBox(self, bbox, wcs, filterName, epoch=None,
515 bboxToSpherePadding=100):
516 """Load reference objects that are within a pixel-based rectangular
517 region.
518
519 This algorithm works by creating a spherical box whose corners
520 correspond to the WCS converted corners of the input bounding box
521 (possibly padded). It then defines a filtering function which looks at
522 the pixel position of the reference objects and accepts only those that
523 lie within the specified bounding box.
524
525 The spherical box region and filtering function are passed to the
526 generic loadRegion method which loads and filters the reference objects
527 from the datastore and returns a single catalog containing the filtered
528 set of reference objects.
529
530 Parameters
531 ----------
532 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D`
533 Box which bounds a region in pixel space.
534 wcs : `lsst.afw.geom.SkyWcs`
535 Wcs object defining the pixel to sky (and inverse) transform for
536 the supplied ``bbox``.
537 filterName : `str`
538 Name of camera filter.
539 epoch : `astropy.time.Time` or `None`, optional
540 Epoch to which to correct proper motion and parallax, or `None`
541 to not apply such corrections.
542 bboxToSpherePadding : `int`, optional
543 Padding to account for translating a set of corners into a
544 spherical (convex) boundary that is certain to encompase the
545 enitre area covered by the box.
546
547 Returns
548 -------
549 output : `lsst.pipe.base.Struct`
550 Results struct with attributes:
551
552 ``refCat``
553 Catalog containing reference objects inside the specified
554 bounding box (padded by self.config.pixelMargin).
555 ``fluxField``
556 Name of the field containing the flux associated with
557 ``filterName``.
558
559 Raises
560 ------
561 RuntimeError
562 Raised if no reference catalogs could be found for the specified
563 region.
564 TypeError
565 Raised if the loaded reference catalogs do not have matching
566 schemas.
567 """
568 paddedBbox = geom.Box2D(bbox)
569 paddedBbox.grow(self.config.pixelMargin)
570 innerSkyRegion, outerSkyRegion, _, _ = self._makeBoxRegion(paddedBbox, wcs, bboxToSpherePadding)
571
572 def _filterFunction(refCat, region):
573 # Perform an initial "pre filter" step based on the refCat coords
574 # and the outerSkyRegion created from the self.config.pixelMargin-
575 # paddedBbox plus an "extra" padding of bboxToSpherePadding and the
576 # raw wcs. This should ensure a large enough projected area on the
577 # sky that accounts for any projection/distortion issues, but small
578 # enough to filter out loaded reference objects that lie well
579 # beyond the projected detector of interest. This step is required
580 # due to the very local nature of the wcs available for the
581 # sky <--> pixel conversions.
582 preFiltFunc = _FilterCatalog(outerSkyRegion)
583 refCat = preFiltFunc(refCat, region)
584
585 # Add columns to the pre-filtered reference catalog relating their
586 # coordinates to equivalent pixel positions for the wcs provided
587 # and use to populate those columns.
588 refCat = self._remapReferenceCatalogSchema(refCat, centroids=True)
589 afwTable.updateRefCentroids(wcs, refCat)
590 # No need to filter the catalog if it is entirely contained in the
591 # region defined by the inner sky region.
592 if innerSkyRegion.contains(region):
593 return refCat
594 # Create a new reference catalog, and populate it only with records
595 # that fall inside the padded bbox.
596 filteredRefCat = type(refCat)(refCat.table)
597 centroidKey = afwTable.Point2DKey(refCat.schema['centroid'])
598 for record in refCat:
599 pixCoords = record[centroidKey]
600 if paddedBbox.contains(geom.Point2D(pixCoords)):
601 filteredRefCat.append(record)
602 return filteredRefCat
603 return self.loadRegion(outerSkyRegion, filterName, filtFunc=_filterFunction, epoch=epoch)
604
605 def loadRegion(self, region, filterName, filtFunc=None, epoch=None):
606 """Load reference objects within a specified region.
607
608 This function loads the DataIds used to construct an instance of this
609 class which intersect or are contained within the specified region. The
610 reference catalogs which intersect but are not fully contained within
611 the input region are further filtered by the specified filter function.
612 This function returns a single source catalog containing all reference
613 objects inside the specified region.
614
615 Parameters
616 ----------
617 region : `lsst.sphgeom.Region`
618 This can be any type that is derived from `lsst.sphgeom.Region` and
619 should define the spatial region for which reference objects are to
620 be loaded.
621 filtFunc : callable or `None`, optional
622 This optional parameter should be a callable object that takes a
623 reference catalog and its corresponding region as parameters,
624 filters the catalog by some criteria and returns the filtered
625 reference catalog. If `None`, an internal filter function is used
626 which filters according to if a reference object falls within the
627 input region.
628 filterName : `str`
629 Name of camera filter.
630 epoch : `astropy.time.Time` or `None`, optional
631 Epoch to which to correct proper motion and parallax, or `None` to
632 not apply such corrections.
633
634 Returns
635 -------
636 output : `lsst.pipe.base.Struct`
637 Results struct with attributes:
638
639 ``refCat``
640 Catalog containing reference objects which intersect the
641 input region, filtered by the specified filter function.
642 ``fluxField``
643 Name of the field containing the flux associated with
644 ``filterName``.
645
646 Raises
647 ------
648 RuntimeError
649 Raised if no reference catalogs could be found for the specified
650 region.
651 TypeError
652 Raised if the loaded reference catalogs do not have matching
653 schemas.
654 """
655 regionLat = region.getBoundingBox().getLat()
656 regionLon = region.getBoundingBox().getLon()
657 self.log.info("Loading reference objects from %s in region bounded by "
658 "[%.8f, %.8f], [%.8f, %.8f] RA Dec",
659 self.name,
660 regionLon.getA().asDegrees(), regionLon.getB().asDegrees(),
661 regionLat.getA().asDegrees(), regionLat.getB().asDegrees())
662 if filtFunc is None:
663 filtFunc = _FilterCatalog(region)
664 # filter out all the regions supplied by the constructor that do not overlap
665 overlapList = []
666 for dataId, refCat in zip(self.dataIds, self.refCats):
667 # SphGeom supports some objects intersecting others, but is not symmetric,
668 # try the intersect operation in both directions
669 try:
670 intersects = dataId.region.intersects(region)
671 except TypeError:
672 intersects = region.intersects(dataId.region)
673
674 if intersects:
675 overlapList.append((dataId, refCat))
676
677 if len(overlapList) == 0:
678 raise RuntimeError("No reference tables could be found for input region")
679
680 firstCat = overlapList[0][1].get()
681 refCat = filtFunc(firstCat, overlapList[0][0].region)
682 trimmedAmount = len(firstCat) - len(refCat)
683
684 # Load in the remaining catalogs
685 for dataId, inputRefCat in overlapList[1:]:
686 tmpCat = inputRefCat.get()
687
688 if tmpCat.schema != firstCat.schema:
689 raise TypeError("Reference catalogs have mismatching schemas")
690
691 filteredCat = filtFunc(tmpCat, dataId.region)
692 refCat.extend(filteredCat)
693 trimmedAmount += len(tmpCat) - len(filteredCat)
694
695 version = getFormatVersionFromRefCat(refCat)
696 if version > LATEST_FORMAT_VERSION:
697 raise ValueError(f"Unsupported refcat format version: {version} > {LATEST_FORMAT_VERSION}.")
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 return pipeBase.Struct(refCat=expandedCat, fluxField=fluxField)
719
720 def loadSkyCircle(self, ctrCoord, radius, filterName, epoch=None):
721 """Load reference objects that lie within a circular region on the sky.
722
723 This method constructs a circular region from an input center and
724 angular radius, loads reference catalogs which are contained in or
725 intersect the circle, and filters reference catalogs which intersect
726 down to objects which lie within the defined circle.
727
728 Parameters
729 ----------
730 ctrCoord : `lsst.geom.SpherePoint`
731 Point defining the center of the circular region.
732 radius : `lsst.geom.Angle`
733 Defines the angular radius of the circular region.
734 filterName : `str`
735 Name of camera filter.
736 epoch : `astropy.time.Time` or `None`, optional
737 Epoch to which to correct proper motion and parallax, or `None` to
738 not apply such corrections.
739
740 Returns
741 -------
742 output : `lsst.pipe.base.Struct`
743 Results struct with attributes:
744
745 ``refCat``
746 Catalog containing reference objects inside the specified
747 search circle.
748 ``fluxField``
749 Name of the field containing the flux associated with
750 ``filterName``.
751 """
752 centerVector = ctrCoord.getVector()
753 sphRadius = sphgeom.Angle(radius.asRadians())
754 circularRegion = sphgeom.Circle(centerVector, sphRadius)
755 return self.loadRegion(circularRegion, filterName, epoch=epoch)
756
757
758def getRefFluxField(schema, filterName):
759 """Get the name of a flux field from a schema.
760
761 Parameters
762 ----------
763 schema : `lsst.afw.table.Schema`
764 Reference catalog schema.
765 filterName : `str`
766 Name of camera filter.
767
768 Returns
769 -------
770 fluxFieldName : `str`
771 Name of flux field.
772
773 Notes
774 -----
775 Return the alias of ``anyFilterMapsToThis``, if present
776 else, return ``*filterName*_camFlux`` if present,
777 else, return ``*filterName*_flux`` if present (camera filter name
778 matches reference filter name), else raise an exception.
779
780 Raises
781 ------
782 RuntimeError
783 Raised if an appropriate field is not found.
784 """
785 if not isinstance(schema, afwTable.Schema):
786 raise RuntimeError("schema=%s is not a schema" % (schema,))
787 try:
788 return schema.getAliasMap().get("anyFilterMapsToThis")
789 except LookupError:
790 pass # try the filterMap next
791
792 fluxFieldList = [filterName + "_camFlux", filterName + "_flux"]
793 for fluxField in fluxFieldList:
794 if fluxField in schema:
795 return fluxField
796
797 raise RuntimeError("Could not find flux field(s) %s" % (", ".join(fluxFieldList)))
798
799
800def getRefFluxKeys(schema, filterName):
801 """Return keys for flux and flux error.
802
803 Parameters
804 ----------
805 schema : `lsst.afw.table.Schema`
806 Reference catalog schema.
807 filterName : `str`
808 Name of camera filter.
809
810 Returns
811 -------
812 keys : `tuple` of (`lsst.afw.table.Key`, `lsst.afw.table.Key`)
813 Two keys:
814
815 - flux key
816 - flux error key, if present, else None
817
818 Raises
819 ------
820 RuntimeError
821 If flux field not found.
822 """
823 fluxField = getRefFluxField(schema, filterName)
824 fluxErrField = fluxField + "Err"
825 fluxKey = schema[fluxField].asKey()
826 try:
827 fluxErrKey = schema[fluxErrField].asKey()
828 except Exception:
829 fluxErrKey = None
830 return (fluxKey, fluxErrKey)
831
832
833def applyProperMotionsImpl(log, catalog, epoch):
834 """Apply proper motion correction to a reference catalog.
835
836 Adjust position and position error in the ``catalog``
837 for proper motion to the specified ``epoch``,
838 modifying the catalog in place.
839
840 Parameters
841 ----------
842 log : `lsst.log.Log` or `logging.getLogger`
843 Log object to write to.
844 catalog : `lsst.afw.table.SimpleCatalog`
845 Catalog of positions, containing:
846
847 - Coordinates, retrieved by the table's coordinate key.
848 - ``coord_raErr`` : Error in Right Ascension (rad).
849 - ``coord_decErr`` : Error in Declination (rad).
850 - ``pm_ra`` : Proper motion in Right Ascension (rad/yr,
851 East positive)
852 - ``pm_raErr`` : Error in ``pm_ra`` (rad/yr), optional.
853 - ``pm_dec`` : Proper motion in Declination (rad/yr,
854 North positive)
855 - ``pm_decErr`` : Error in ``pm_dec`` (rad/yr), optional.
856 - ``epoch`` : Mean epoch of object (an astropy.time.Time)
857 epoch : `astropy.time.Time`
858 Epoch to which to correct proper motion.
859 """
860 if "epoch" not in catalog.schema or "pm_ra" not in catalog.schema or "pm_dec" not in catalog.schema:
861 log.warning("Proper motion correction not available from catalog")
862 return
863 if not catalog.isContiguous():
864 raise RuntimeError("Catalog must be contiguous")
865 catEpoch = astropy.time.Time(catalog["epoch"], scale="tai", format="mjd")
866 log.info("Correcting reference catalog for proper motion to %r", epoch)
867 # Use `epoch.tai` to make sure the time difference is in TAI
868 timeDiffsYears = (epoch.tai - catEpoch).to(astropy.units.yr).value
869 coordKey = catalog.table.getCoordKey()
870 # Compute the offset of each object due to proper motion
871 # as components of the arc of a great circle along RA and Dec
872 pmRaRad = catalog["pm_ra"]
873 pmDecRad = catalog["pm_dec"]
874 offsetsRaRad = pmRaRad*timeDiffsYears
875 offsetsDecRad = pmDecRad*timeDiffsYears
876 # Compute the corresponding bearing and arc length of each offset
877 # due to proper motion, and apply the offset.
878 # The factor of 1e6 for computing bearing is intended as
879 # a reasonable scale for typical values of proper motion
880 # in order to avoid large errors for small values of proper motion;
881 # using the offsets is another option, but it can give
882 # needlessly large errors for short duration.
883 offsetBearingsRad = numpy.arctan2(offsetsDecRad*1e6, offsetsRaRad*1e6)
884 offsetAmountsRad = numpy.hypot(offsetsRaRad, offsetsDecRad)
885 for record, bearingRad, amountRad in zip(catalog, offsetBearingsRad, offsetAmountsRad):
886 record.set(coordKey,
887 record.get(coordKey).offset(bearing=bearingRad*geom.radians,
888 amount=amountRad*geom.radians))
889 # TODO DM-36979: this needs to incorporate the full covariance!
890 # Increase error in RA and Dec based on error in proper motion
891 if "coord_raErr" in catalog.schema:
892 catalog["coord_raErr"] = numpy.hypot(catalog["coord_raErr"],
893 catalog["pm_raErr"]*timeDiffsYears)
894 if "coord_decErr" in catalog.schema:
895 catalog["coord_decErr"] = numpy.hypot(catalog["coord_decErr"],
896 catalog["pm_decErr"]*timeDiffsYears)
int max
Tag types used to declare specialized field types.
Definition misc.h:31
Defines the fields and offsets for a table.
Definition Schema.h:51
A mapping between the keys of two Schemas, used to copy data between them.
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
Class for storing ordered metadata with comments.
A floating-point coordinate rectangle geometry.
Definition Box.h:413
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
loadRegion(self, region, filterName, filtFunc=None, epoch=None)
_addFluxAliases(schema, anyFilterMapsToThis=None, filterMap=None)
getMetadataBox(self, bbox, wcs, filterName, epoch=None, bboxToSpherePadding=100)
__init__(self, dataIds, refCats, name=None, log=None, config=None)
loadPixelBox(self, bbox, wcs, filterName, epoch=None, bboxToSpherePadding=100)
_remapReferenceCatalogSchema(refCat, *anyFilterMapsToThis=None, filterMap=None, centroids=False)
loadSkyCircle(self, ctrCoord, radius, filterName, epoch=None)
Angle represents an angle in radians.
Definition Angle.h:50
Circle is a circular region on the unit sphere that contains its boundary.
Definition Circle.h:53
ConvexPolygon is a closed convex polygon on the unit sphere.
daf::base::PropertySet * set
Definition fits.cc:931
void updateRefCentroids(geom::SkyWcs const &wcs, ReferenceCollection &refList)
Update centroids in a collection of reference objects.
Definition wcsUtils.cc:73