Loading [MathJax]/extensions/tex2jax.js
LSST Applications 28.0.0,g1653933729+a8ce1bb630,g1a997c3884+a8ce1bb630,g28da252d5a+5bd70b7e6d,g2bbee38e9b+638fca75ac,g2bc492864f+638fca75ac,g3156d2b45e+07302053f8,g347aa1857d+638fca75ac,g35bb328faa+a8ce1bb630,g3a166c0a6a+638fca75ac,g3e281a1b8c+7bbb0b2507,g4005a62e65+17cd334064,g414038480c+5b5cd4fff3,g41af890bb2+4ffae9de63,g4e1a3235cc+0f1912dca3,g6249c6f860+3c3976f90c,g80478fca09+46aba80bd6,g82479be7b0+77990446f6,g858d7b2824+78ba4d1ce1,g89c8672015+f667a5183b,g9125e01d80+a8ce1bb630,ga5288a1d22+2a6264e9ca,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gc22bb204ba+78ba4d1ce1,gc28159a63d+638fca75ac,gcf0d15dbbd+32ddb6096f,gd6b7c0dfd1+3e339405e9,gda3e153d99+78ba4d1ce1,gda6a2b7d83+32ddb6096f,gdaeeff99f8+1711a396fd,gdd5a9049c5+b18c39e5e3,ge2409df99d+a5e4577cdc,ge33fd446bb+78ba4d1ce1,ge79ae78c31+638fca75ac,gf0baf85859+64e8883e75,gf5289d68f6+e1b046a8d7,gfa443fc69c+91d9ed1ecf,gfda6b12a05+8419469a56
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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 self.log.debug("Trimmed %d refCat objects lying outside padded region, leaving %d",
696 trimmedAmount, len(refCat))
697 self.log.info("Loaded %d reference objects", len(refCat))
698
699 # Ensure that the loaded reference catalog is continuous in memory
700 if not refCat.isContiguous():
701 refCat = refCat.copy(deep=True)
702
703 self.applyProperMotions(refCat, epoch)
704
705 expandedCat = self._remapReferenceCatalogSchema(refCat,
706 anyFilterMapsToThis=self.config.anyFilterMapsToThis,
707 filterMap=self.config.filterMap)
708
709 # Ensure that the returned reference catalog is continuous in memory
710 if not expandedCat.isContiguous():
711 expandedCat = expandedCat.copy(deep=True)
712
713 fluxField = getRefFluxField(expandedCat.schema, filterName)
714
715 if expandedCat.schema[fluxField].asField().getUnits() != "nJy":
716 # if the flux field is not in nJy, check the refcat format version
717 version = getFormatVersionFromRefCat(refCat)
718 if version > LATEST_FORMAT_VERSION:
719 raise ValueError(f"Unsupported refcat format version: {version} > {LATEST_FORMAT_VERSION}.")
720
721 return pipeBase.Struct(refCat=expandedCat, fluxField=fluxField)
722
723 def loadSkyCircle(self, ctrCoord, radius, filterName, epoch=None):
724 """Load reference objects that lie within a circular region on the sky.
725
726 This method constructs a circular region from an input center and
727 angular radius, loads reference catalogs which are contained in or
728 intersect the circle, and filters reference catalogs which intersect
729 down to objects which lie within the defined circle.
730
731 Parameters
732 ----------
733 ctrCoord : `lsst.geom.SpherePoint`
734 Point defining the center of the circular region.
735 radius : `lsst.geom.Angle`
736 Defines the angular radius of the circular region.
737 filterName : `str`
738 Name of camera filter.
739 epoch : `astropy.time.Time` or `None`, optional
740 Epoch to which to correct proper motion and parallax, or `None` to
741 not apply such corrections.
742
743 Returns
744 -------
745 output : `lsst.pipe.base.Struct`
746 Results struct with attributes:
747
748 ``refCat``
749 Catalog containing reference objects inside the specified
750 search circle.
751 ``fluxField``
752 Name of the field containing the flux associated with
753 ``filterName``.
754 """
755 centerVector = ctrCoord.getVector()
756 sphRadius = sphgeom.Angle(radius.asRadians())
757 circularRegion = sphgeom.Circle(centerVector, sphRadius)
758 return self.loadRegion(circularRegion, filterName, epoch=epoch)
759
760
761def getRefFluxField(schema, filterName):
762 """Get the name of a flux field from a schema.
763
764 Parameters
765 ----------
766 schema : `lsst.afw.table.Schema`
767 Reference catalog schema.
768 filterName : `str`
769 Name of camera filter.
770
771 Returns
772 -------
773 fluxFieldName : `str`
774 Name of flux field.
775
776 Notes
777 -----
778 Return the alias of ``anyFilterMapsToThis``, if present
779 else, return ``*filterName*_camFlux`` if present,
780 else, return ``*filterName*_flux`` if present (camera filter name
781 matches reference filter name), else raise an exception.
782
783 Raises
784 ------
785 RuntimeError
786 Raised if an appropriate field is not found.
787 """
788 if not isinstance(schema, afwTable.Schema):
789 raise RuntimeError("schema=%s is not a schema" % (schema,))
790 try:
791 return schema.getAliasMap().get("anyFilterMapsToThis")
792 except LookupError:
793 pass # try the filterMap next
794
795 fluxFieldList = [filterName + "_camFlux", filterName + "_flux"]
796 for fluxField in fluxFieldList:
797 if fluxField in schema:
798 return fluxField
799
800 raise RuntimeError("Could not find flux field(s) %s" % (", ".join(fluxFieldList)))
801
802
803def getRefFluxKeys(schema, filterName):
804 """Return keys for flux and flux error.
805
806 Parameters
807 ----------
808 schema : `lsst.afw.table.Schema`
809 Reference catalog schema.
810 filterName : `str`
811 Name of camera filter.
812
813 Returns
814 -------
815 keys : `tuple` of (`lsst.afw.table.Key`, `lsst.afw.table.Key`)
816 Two keys:
817
818 - flux key
819 - flux error key, if present, else None
820
821 Raises
822 ------
823 RuntimeError
824 If flux field not found.
825 """
826 fluxField = getRefFluxField(schema, filterName)
827 fluxErrField = fluxField + "Err"
828 fluxKey = schema[fluxField].asKey()
829 try:
830 fluxErrKey = schema[fluxErrField].asKey()
831 except Exception:
832 fluxErrKey = None
833 return (fluxKey, fluxErrKey)
834
835
836def applyProperMotionsImpl(log, catalog, epoch):
837 """Apply proper motion correction to a reference catalog.
838
839 Adjust position and position error in the ``catalog``
840 for proper motion to the specified ``epoch``,
841 modifying the catalog in place.
842
843 Parameters
844 ----------
845 log : `lsst.log.Log` or `logging.Logger`
846 Log object to write to.
847 catalog : `lsst.afw.table.SimpleCatalog`
848 Catalog of positions, containing:
849
850 - Coordinates, retrieved by the table's coordinate key.
851 - ``coord_raErr`` : Error in Right Ascension (rad).
852 - ``coord_decErr`` : Error in Declination (rad).
853 - ``pm_ra`` : Proper motion in Right Ascension (rad/yr,
854 East positive)
855 - ``pm_raErr`` : Error in ``pm_ra`` (rad/yr), optional.
856 - ``pm_dec`` : Proper motion in Declination (rad/yr,
857 North positive)
858 - ``pm_decErr`` : Error in ``pm_dec`` (rad/yr), optional.
859 - ``epoch`` : Mean epoch of object (an astropy.time.Time)
860 epoch : `astropy.time.Time`
861 Epoch to which to correct proper motion.
862 """
863 if "epoch" not in catalog.schema or "pm_ra" not in catalog.schema or "pm_dec" not in catalog.schema:
864 log.warning("Proper motion correction not available from catalog")
865 return
866 if not catalog.isContiguous():
867 raise RuntimeError("Catalog must be contiguous")
868 catEpoch = astropy.time.Time(catalog["epoch"], scale="tai", format="mjd")
869 log.info("Correcting reference catalog for proper motion to %r", epoch)
870 # Use `epoch.tai` to make sure the time difference is in TAI
871 timeDiffsYears = (epoch.tai - catEpoch).to(astropy.units.yr).value
872 coordKey = catalog.table.getCoordKey()
873 # Compute the offset of each object due to proper motion
874 # as components of the arc of a great circle along RA and Dec
875 pmRaRad = catalog["pm_ra"]
876 pmDecRad = catalog["pm_dec"]
877 offsetsRaRad = pmRaRad*timeDiffsYears
878 offsetsDecRad = pmDecRad*timeDiffsYears
879 # Compute the corresponding bearing and arc length of each offset
880 # due to proper motion, and apply the offset.
881 # The factor of 1e6 for computing bearing is intended as
882 # a reasonable scale for typical values of proper motion
883 # in order to avoid large errors for small values of proper motion;
884 # using the offsets is another option, but it can give
885 # needlessly large errors for short duration.
886 offsetBearingsRad = numpy.arctan2(offsetsDecRad*1e6, offsetsRaRad*1e6)
887 offsetAmountsRad = numpy.hypot(offsetsRaRad, offsetsDecRad)
888 for record, bearingRad, amountRad in zip(catalog, offsetBearingsRad, offsetAmountsRad):
889 record.set(coordKey,
890 record.get(coordKey).offset(bearing=bearingRad*geom.radians,
891 amount=amountRad*geom.radians))
892 # TODO DM-36979: this needs to incorporate the full covariance!
893 # Increase error in RA and Dec based on error in proper motion
894 if "coord_raErr" in catalog.schema:
895 catalog["coord_raErr"] = numpy.hypot(catalog["coord_raErr"],
896 catalog["pm_raErr"]*timeDiffsYears)
897 if "coord_decErr" in catalog.schema:
898 catalog["coord_decErr"] = numpy.hypot(catalog["coord_decErr"],
899 catalog["pm_decErr"]*timeDiffsYears)
int max
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)
loadSkyCircle(self, ctrCoord, radius, filterName, epoch=None)
_remapReferenceCatalogSchema(refCat, *, anyFilterMapsToThis=None, filterMap=None, centroids=False)
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:54
ConvexPolygon is a closed convex polygon on the unit sphere.
void updateRefCentroids(geom::SkyWcs const &wcs, ReferenceCollection &refList)
Update centroids in a collection of reference objects.
Definition wcsUtils.cc:73