LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
loadReferenceCatalog.py
Go to the documentation of this file.
1 # See COPYRIGHT file at the top of the source tree.
2 #
3 # This file is part of pipe_tasks.
4 #
5 # Developed for the LSST Data Management System.
6 # This product includes software developed by the LSST Project
7 # (https://www.lsst.org).
8 # See the COPYRIGHT file at the top-level directory of this distribution
9 # for details of code ownership.
10 #
11 # This program is free software: you can redistribute it and/or modify
12 # it under the terms of the GNU General Public License as published by
13 # the Free Software Foundation, either version 3 of the License, or
14 # (at your option) any later version.
15 #
16 # This program is distributed in the hope that it will be useful,
17 # but WITHOUT ANY WARRANTY; without even the implied warranty of
18 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 # GNU General Public License for more details.
20 #
21 # You should have received a copy of the GNU General Public License
22 # along with this program. If not, see <https://www.gnu.org/licenses/>.
23 """Load a full reference catalog in numpy/table/dataframe format.
24 
25 This task will load multi-band reference objects, apply a reference selector,
26 and apply color terms.
27 """
28 import numpy as np
29 from astropy import units
30 
31 import lsst.pex.config as pexConfig
32 import lsst.pipe.base as pipeBase
33 from lsst.meas.algorithms import LoadIndexedReferenceObjectsTask, ReferenceSourceSelectorTask
34 from lsst.meas.algorithms import getRefFluxField
35 from lsst.pipe.tasks.colorterms import ColortermLibrary
36 from lsst.afw.image import abMagErrFromFluxErr
37 from lsst.meas.algorithms import ReferenceObjectLoader, LoadReferenceObjectsConfig
38 
39 import lsst.geom
40 
41 __all__ = ['LoadReferenceCatalogConfig', 'LoadReferenceCatalogTask']
42 
43 
44 class LoadReferenceCatalogConfig(pexConfig.Config):
45  """Config for LoadReferenceCatalogTask"""
46  refObjLoader = pexConfig.ConfigurableField(
47  target=LoadIndexedReferenceObjectsTask,
48  doc="Reference object loader for photometry",
49  )
50  doApplyColorTerms = pexConfig.Field(
51  doc=("Apply photometric color terms to reference stars? "
52  "Requires that colorterms be set to a ColorTermLibrary"),
53  dtype=bool,
54  default=True
55  )
56  colorterms = pexConfig.ConfigField(
57  doc="Library of photometric reference catalog name to color term dict.",
58  dtype=ColortermLibrary,
59  )
60  referenceSelector = pexConfig.ConfigurableField(
61  target=ReferenceSourceSelectorTask,
62  doc="Selection of reference sources",
63  )
64  doReferenceSelection = pexConfig.Field(
65  doc="Run the reference selector on the reference catalog?",
66  dtype=bool,
67  default=True
68  )
69 
70  def validate(self):
71  super().validate()
72  if self.doApplyColorTermsdoApplyColorTerms and len(self.colortermscolorterms.data) == 0:
73  msg = "applyColorTerms=True requires the `colorterms` field be set to a ColortermLibrary."
74  raise pexConfig.FieldValidationError(LoadReferenceCatalogConfig.colorterms, self, msg)
75 
76 
77 class LoadReferenceCatalogTask(pipeBase.Task):
78  """Load multi-band reference objects from a reference catalog.
79 
80  Parameters
81  ----------
82  dataIds : iterable of `lsst.daf.butler.dataId`, optional
83  An iterable object of dataIds which point to reference catalogs
84  in a Gen3 repository. Required for Gen3.
85  refCats : iterable of `lsst.daf.butler.DeferredDatasetHandle`, optional
86  An iterable object of dataset refs for reference catalogs in
87  a Gen3 repository. Required for Gen3.
88  butler : `lsst.daf.persistence.Butler`, optional
89  A Gen2 butler. Required for Gen2.
90  """
91  ConfigClass = LoadReferenceCatalogConfig
92  _DefaultName = "loadReferenceCatalog"
93 
94  def __init__(self, dataIds=None, refCats=None, butler=None, **kwargs):
95  if dataIds is not None and refCats is not None:
96  pipeBase.Task.__init__(self, **kwargs)
97  refConfig = self.config.refObjLoader
98  self.refObjLoaderrefObjLoader = ReferenceObjectLoader(dataIds=dataIds,
99  refCats=refCats,
100  config=refConfig,
101  log=self.log)
102  elif butler is not None:
103  raise RuntimeError("LoadReferenceCatalogTask does not support Gen2 Butlers.")
104  else:
105  raise RuntimeError("Must instantiate LoadReferenceCatalogTask with "
106  "dataIds and refCats (Gen3)")
107 
108  if self.config.doReferenceSelection:
109  self.makeSubtask('referenceSelector')
110  self._fluxFilters_fluxFilters = None
111  self._fluxFields_fluxFields = None
112  self._referenceFilter_referenceFilter = None
113 
114  def getPixelBoxCatalog(self, bbox, wcs, filterList, epoch=None,
115  bboxToSpherePadding=None):
116  """Get a multi-band reference catalog by specifying a bounding box and WCS.
117 
118  The catalog will be in `numpy.ndarray`, with positions proper-motion
119  corrected to "epoch" (if specified, and if the reference catalog has
120  proper motions); sources cut on a reference selector (if
121  "config.doReferenceSelection = True"); and color-terms applied (if
122  "config.doApplyColorTerms = True").
123 
124  The format of the reference catalog will be of the format:
125 
126  dtype = [('ra', 'np.float64'),
127  ('dec', 'np.float64'),
128  ('refMag', 'np.float32', (len(filterList), )),
129  ('refMagErr', 'np.float32', (len(filterList), ))]
130 
131  Reference magnitudes (AB) and errors will be 99 for non-detections
132  in a given band.
133 
134  Parameters
135  ----------
136  bbox : `lsst.geom.Box2I`
137  Box which bounds a region in pixel space.
138  wcs : `lsst.afw.geom.SkyWcs`
139  Wcs object defining the pixel to sky (and reverse) transform for
140  the supplied bbox.
141  filterList : `List` [ `str` ]
142  List of camera physicalFilter names to retrieve magnitudes.
143  epoch : `astropy.time.Time`, optional
144  Epoch to which to correct proper motion and parallax
145  (if available), or `None` to not apply such corrections.
146  bboxToSpherePadding : `int`, optional
147  Padding to account for translating a set of corners into a
148  spherical (convex) boundary that is certain to encompass the
149  entire area covered by the bbox.
150 
151  Returns
152  -------
153  refCat : `numpy.ndarray`
154  Reference catalog.
155  """
156  # Check if we have previously cached values for the fluxFields
157  if self._fluxFilters_fluxFilters is None or self._fluxFilters_fluxFilters != filterList:
158  center = wcs.pixelToSky(bbox.getCenter())
159  self._determineFluxFields_determineFluxFields(center, filterList)
160 
161  skyBox = self.refObjLoaderrefObjLoader.loadPixelBox(bbox, wcs, self._referenceFilter_referenceFilter,
162  epoch=epoch,
163  bboxToSpherePadding=bboxToSpherePadding)
164 
165  if not skyBox.refCat.isContiguous():
166  refCat = skyBox.refCat.copy(deep=True)
167  else:
168  refCat = skyBox.refCat
169 
170  return self._formatCatalog_formatCatalog(refCat, filterList)
171 
172  def getSkyCircleCatalog(self, center, radius, filterList, epoch=None,
173  catalogFormat='numpy'):
174  """Get a multi-band reference catalog by specifying a center and radius.
175 
176  The catalog will be in `numpy.ndarray`, with positions proper-motion
177  corrected to "epoch" (if specified, and if the reference catalog has
178  proper motions); sources cut on a reference selector (if
179  "config.doReferenceSelection = True"); and color-terms applied (if
180  "config.doApplyColorTerms = True").
181 
182  The format of the reference catalog will be of the format:
183 
184  dtype = [('ra', 'np.float64'),
185  ('dec', 'np.float64'),
186  ('refMag', 'np.float32', (len(filterList), )),
187  ('refMagErr', 'np.float32', (len(filterList), ))]
188 
189  Reference magnitudes (AB) and errors will be 99 for non-detections
190  in a given band.
191 
192  Parameters
193  ----------
194  center : `lsst.geom.SpherePoint`
195  Point defining the center of the circular region.
196  radius : `lsst.geom.Angle`
197  Defines the angular radius of the circular region.
198  filterList : `List` [ `str` ]
199  List of camera physicalFilter names to retrieve magnitudes.
200  epoch : `astropy.time.Time`, optional
201  Epoch to which to correct proper motion and parallax
202  (if available), or `None` to not apply such corrections.
203 
204  Parameters
205  ----------
206  refCat : `numpy.ndarray`
207  Reference catalog.
208  """
209  # Check if we have previously cached values for the fluxFields
210  if self._fluxFilters_fluxFilters is None or self._fluxFilters_fluxFilters != filterList:
211  self._determineFluxFields_determineFluxFields(center, filterList)
212 
213  skyCircle = self.refObjLoaderrefObjLoader.loadSkyCircle(center, radius,
214  self._referenceFilter_referenceFilter,
215  epoch=epoch)
216 
217  if not skyCircle.refCat.isContiguous():
218  refCat = skyCircle.refCat.copy(deep=True)
219  else:
220  refCat = skyCircle.refCat
221 
222  return self._formatCatalog_formatCatalog(refCat, filterList)
223 
224  def _formatCatalog(self, refCat, filterList):
225  """Format a reference afw table into the final format.
226 
227  This method applies reference selections and color terms as specified
228  by the config.
229 
230  Parameters
231  ----------
232  refCat : `lsst.afw.table.SourceCatalog`
233  Reference catalog in afw format.
234  filterList : `list` [`str`]
235  List of camera physicalFilter names to apply color terms.
236 
237  Returns
238  -------
239  refCat : `numpy.ndarray`
240  Reference catalog.
241  """
242  if self.config.doReferenceSelection:
243  goodSources = self.referenceSelector.selectSources(refCat)
244  selected = goodSources.selected
245  else:
246  selected = np.ones(len(refCat), dtype=bool)
247 
248  npRefCat = np.zeros(np.sum(selected), dtype=[('ra', 'f8'),
249  ('dec', 'f8'),
250  ('refMag', 'f4', (len(filterList), )),
251  ('refMagErr', 'f4', (len(filterList), ))])
252 
253  if npRefCat.size == 0:
254  # Return an empty catalog if we don't have any selected sources.
255  return npRefCat
256 
257  # Natively "coord_ra" and "coord_dec" are stored in radians.
258  # Doing this as an array rather than by row with the coord access is
259  # approximately 600x faster.
260  npRefCat['ra'] = np.rad2deg(refCat['coord_ra'][selected])
261  npRefCat['dec'] = np.rad2deg(refCat['coord_dec'][selected])
262 
263  # Default (unset) values are 99.0
264  npRefCat['refMag'][:, :] = 99.0
265  npRefCat['refMagErr'][:, :] = 99.0
266 
267  if self.config.doApplyColorTerms:
268  if isinstance(self.refObjLoaderrefObjLoader, ReferenceObjectLoader):
269  # Gen3
270  refCatName = self.refObjLoaderrefObjLoader.config.value.ref_dataset_name
271  else:
272  # Gen2
273  refCatName = self.refObjLoaderrefObjLoader.ref_dataset_name
274 
275  for i, (filterName, fluxField) in enumerate(zip(self._fluxFilters_fluxFilters, self._fluxFields_fluxFields)):
276  if fluxField is None:
277  # There is no matching reference band.
278  # This will leave the column filled with 99s
279  continue
280  self.log.debug("Applying color terms for filterName='%s'", filterName)
281 
282  colorterm = self.config.colorterms.getColorterm(filterName, refCatName, doRaise=True)
283 
284  refMag, refMagErr = colorterm.getCorrectedMagnitudes(refCat)
285 
286  # nan_to_num replaces nans with zeros, and this ensures
287  # that we select magnitudes that both filter out nans and are
288  # not very large (corresponding to very small fluxes), as "99"
289  # is a commen sentinel for illegal magnitudes.
290  good, = np.where((np.nan_to_num(refMag[selected], nan=99.0) < 90.0)
291  & (np.nan_to_num(refMagErr[selected], nan=99.0) < 90.0)
292  & (np.nan_to_num(refMagErr[selected]) > 0.0))
293 
294  npRefCat['refMag'][good, i] = refMag[selected][good]
295  npRefCat['refMagErr'][good, i] = refMagErr[selected][good]
296  else:
297  # No color terms to apply
298  for i, (filterName, fluxField) in enumerate(zip(self._fluxFilters_fluxFilters, self._fluxFields_fluxFields)):
299  # nan_to_num replaces nans with zeros, and this ensures that
300  # we select fluxes that both filter out nans and are positive.
301  good, = np.where((np.nan_to_num(refCat[fluxField][selected]) > 0.0)
302  & (np.nan_to_num(refCat[fluxField+'Err'][selected]) > 0.0))
303  refMag = (refCat[fluxField][selected][good]*units.nJy).to_value(units.ABmag)
304  refMagErr = abMagErrFromFluxErr(refCat[fluxField+'Err'][selected][good],
305  refCat[fluxField][selected][good])
306  npRefCat['refMag'][good, i] = refMag
307  npRefCat['refMagErr'][good, i] = refMagErr
308 
309  return npRefCat
310 
311  def _determineFluxFields(self, center, filterList):
312  """Determine the flux field names for a reference catalog.
313 
314  This method sets self._fluxFields, self._referenceFilter.
315 
316  Parameters
317  ----------
318  center : `lsst.geom.SpherePoint`
319  The center around which to load test sources.
320  filterList : `list` [`str`]
321  List of camera physicalFilter names.
322  """
323  # Search for a good filter to use to load the reference catalog
324  # via the refObjLoader task which requires a valid filterName
325  foundReferenceFilter = False
326 
327  # Store the original config
328  _config = self.refObjLoaderrefObjLoader.config
329 
330  configTemp = LoadReferenceObjectsConfig()
331  configIntersection = {k: getattr(self.refObjLoaderrefObjLoader.config, k)
332  for k, v in self.refObjLoaderrefObjLoader.config.toDict().items()
333  if (k in configTemp.keys() and k != "connections")}
334  # We must turn off the proper motion checking to find the refFilter.
335  configIntersection['requireProperMotion'] = False
336  configTemp.update(**configIntersection)
337 
338  self.refObjLoaderrefObjLoader.config = configTemp
339 
340  for filterName in filterList:
341  if self.config.refObjLoader.anyFilterMapsToThis is not None:
342  refFilterName = self.config.refObjLoader.anyFilterMapsToThis
343  else:
344  refFilterName = self.config.refObjLoader.filterMap.get(filterName)
345  if refFilterName is None:
346  continue
347  try:
348  results = self.refObjLoaderrefObjLoader.loadSkyCircle(center,
349  0.05*lsst.geom.degrees,
350  refFilterName)
351  foundReferenceFilter = True
352  self._referenceFilter_referenceFilter = refFilterName
353  break
354  except RuntimeError as err:
355  # This just means that the filterName wasn't listed
356  # in the reference catalog. This is okay.
357  if 'not find flux' in err.args[0]:
358  # The filterName wasn't listed in the reference catalog.
359  # This is not a fatal failure (yet)
360  pass
361  else:
362  raise err
363 
364  self.refObjLoaderrefObjLoader.config = _config
365 
366  if not foundReferenceFilter:
367  raise RuntimeError("Could not find any valid flux field(s) %s" %
368  (", ".join(filterList)))
369 
370  # Record self._fluxFilters for checks on subsequent calls
371  self._fluxFilters_fluxFilters = filterList
372 
373  # Retrieve all the fluxField names
374  self._fluxFields_fluxFields = []
375  for filterName in filterList:
376  fluxField = None
377 
378  if self.config.refObjLoader.anyFilterMapsToThis is not None:
379  refFilterName = self.config.refObjLoader.anyFilterMapsToThis
380  else:
381  refFilterName = self.config.refObjLoader.filterMap.get(filterName)
382 
383  if refFilterName is not None:
384  try:
385  fluxField = getRefFluxField(results.refCat.schema, filterName=refFilterName)
386  except RuntimeError:
387  # This flux field isn't available. Set to None
388  fluxField = None
389 
390  if fluxField is None:
391  self.log.warning('No reference flux field for camera filter %s', filterName)
392 
393  self._fluxFields_fluxFields.append(fluxField)
std::vector< SchemaItem< Flag > > * items
def getPixelBoxCatalog(self, bbox, wcs, filterList, epoch=None, bboxToSpherePadding=None)
def __init__(self, dataIds=None, refCats=None, butler=None, **kwargs)
def getSkyCircleCatalog(self, center, radius, filterList, epoch=None, catalogFormat='numpy')
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
double abMagErrFromFluxErr(double fluxErr, double flux)
Compute AB magnitude error from flux and flux error in Janskys.
Definition: Calib.h:52
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.