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
astrometry.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2016 AURA/LSST.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <https://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 __all__ = ["AstrometryConfig", "AstrometryTask"]
24 
25 import numpy as np
26 from astropy import units
27 import scipy.stats
28 
29 import lsst.pex.config as pexConfig
30 import lsst.pipe.base as pipeBase
31 from lsst.utils.timer import timeMethod
32 from .ref_match import RefMatchTask, RefMatchConfig
33 from .fitTanSipWcs import FitTanSipWcsTask
34 from .display import displayAstrometry
35 
36 
38  """Config for AstrometryTask.
39  """
40  wcsFitter = pexConfig.ConfigurableField(
41  target=FitTanSipWcsTask,
42  doc="WCS fitter",
43  )
44  forceKnownWcs = pexConfig.Field(
45  dtype=bool,
46  doc="If True then load reference objects and match sources but do not fit a WCS; "
47  "this simply controls whether 'run' calls 'solve' or 'loadAndMatch'",
48  default=False,
49  )
50  maxIter = pexConfig.RangeField(
51  doc="maximum number of iterations of match sources and fit WCS"
52  "ignored if not fitting a WCS",
53  dtype=int,
54  default=3,
55  min=1,
56  )
57  minMatchDistanceArcSec = pexConfig.RangeField(
58  doc="the match distance below which further iteration is pointless (arcsec); "
59  "ignored if not fitting a WCS",
60  dtype=float,
61  default=0.001,
62  min=0,
63  )
64  maxMeanDistanceArcsec = pexConfig.RangeField(
65  doc="Maximum mean on-sky distance (in arcsec) between matched source and rerference "
66  "objects post-fit. A mean distance greater than this threshold raises a TaskError "
67  "and the WCS fit is considered a failure. The default is set to the maximum tolerated "
68  "by the external global calibration (e.g. jointcal) step for conceivable recovery. "
69  "Appropriate value will be dataset and workflow dependent.",
70  dtype=float,
71  default=0.5,
72  min=0,
73  )
74  doMagnitudeOutlierRejection = pexConfig.Field(
75  dtype=bool,
76  doc=("If True then a rough zeropoint will be computed from matched sources "
77  "and outliers will be rejected in the iterations."),
78  default=False,
79  )
80  magnitudeOutlierRejectionNSigma = pexConfig.Field(
81  dtype=float,
82  doc=("Number of sigma (measured from the distribution) in magnitude "
83  "for a potential reference/source match to be rejected during "
84  "iteration."),
85  default=3.0,
86  )
87 
88  def setDefaults(self):
89  # Override the default source selector for astrometry tasks
90  self.sourceFluxTypesourceFluxTypesourceFluxType = "Ap"
91 
92  self.sourceSelectorsourceSelector.name = "matcher"
93  self.sourceSelectorsourceSelector["matcher"].sourceFluxType = self.sourceFluxTypesourceFluxTypesourceFluxType
94 
95  # Note that if the matcher is MatchOptimisticBTask, then the
96  # default should be self.sourceSelector['matcher'].excludePixelFlags = False
97  # However, there is no way to do this automatically.
98 
99 
101  """Match an input source catalog with objects from a reference catalog and
102  solve for the WCS.
103 
104  This task is broken into two main subasks: matching and WCS fitting which
105  are very interactive. The matching here can be considered in part a first
106  pass WCS fitter due to the fitter's sensitivity to outliers.
107 
108  Parameters
109  ----------
110  refObjLoader : `lsst.meas.algorithms.ReferenceLoader`
111  A reference object loader object
112  schema : `lsst.afw.table.Schema`
113  Used to set "calib_astrometry_used" flag in output source catalog.
114  **kwargs
115  additional keyword arguments for pipe_base
116  `lsst.pipe.base.Task.__init__`
117  """
118  ConfigClass = AstrometryConfig
119  _DefaultName = "astrometricSolver"
120 
121  def __init__(self, refObjLoader, schema=None, **kwargs):
122  RefMatchTask.__init__(self, refObjLoader, **kwargs)
123 
124  if schema is not None:
125  self.usedKeyusedKey = schema.addField("calib_astrometry_used", type="Flag",
126  doc="set if source was used in astrometric calibration")
127  else:
128  self.usedKeyusedKey = None
129 
130  self.makeSubtask("wcsFitter")
131 
132  @timeMethod
133  def run(self, sourceCat, exposure):
134  """Load reference objects, match sources and optionally fit a WCS.
135 
136  This is a thin layer around solve or loadAndMatch, depending on
137  config.forceKnownWcs.
138 
139  Parameters
140  ----------
141  exposure : `lsst.afw.image.Exposure`
142  exposure whose WCS is to be fit
143  The following are read only:
144 
145  - bbox
146  - photoCalib (may be absent)
147  - filter (may be unset)
148  - detector (if wcs is pure tangent; may be absent)
149 
150  The following are updated:
151 
152  - wcs (the initial value is used as an initial guess, and is
153  required)
154 
155  sourceCat : `lsst.afw.table.SourceCatalog`
156  catalog of sources detected on the exposure
157 
158  Returns
159  -------
160  result : `lsst.pipe.base.Struct`
161  with these fields:
162 
163  - ``refCat`` : reference object catalog of objects that overlap the
164  exposure (with some margin) (`lsst.afw.table.SimpleCatalog`).
165  - ``matches`` : astrometric matches
166  (`list` of `lsst.afw.table.ReferenceMatch`).
167  - ``scatterOnSky`` : median on-sky separation between reference
168  objects and sources in "matches"
169  (`lsst.afw.geom.Angle`) or `None` if config.forceKnownWcs True
170  - ``matchMeta`` : metadata needed to unpersist matches
171  (`lsst.daf.base.PropertyList`)
172  """
173  if self.refObjLoaderrefObjLoader is None:
174  raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
175  if self.config.forceKnownWcs:
176  res = self.loadAndMatchloadAndMatch(exposure=exposure, sourceCat=sourceCat)
177  res.scatterOnSky = None
178  else:
179  res = self.solvesolve(exposure=exposure, sourceCat=sourceCat)
180  return res
181 
182  @timeMethod
183  def solve(self, exposure, sourceCat):
184  """Load reference objects overlapping an exposure, match to sources and
185  fit a WCS
186 
187  Returns
188  -------
189  result : `lsst.pipe.base.Struct`
190  Result struct with components:
191 
192  - ``refCat`` : reference object catalog of objects that overlap the
193  exposure (with some margin) (`lsst::afw::table::SimpleCatalog`).
194  - ``matches`` : astrometric matches
195  (`list` of `lsst.afw.table.ReferenceMatch`).
196  - ``scatterOnSky`` : median on-sky separation between reference
197  objects and sources in "matches" (`lsst.geom.Angle`)
198  - ``matchMeta`` : metadata needed to unpersist matches
199  (`lsst.daf.base.PropertyList`)
200 
201  Raises
202  ------
203  TaskError
204  If the measured mean on-sky distance between the matched source and
205  reference objects is greater than
206  ``self.config.maxMeanDistanceArcsec``.
207 
208  Notes
209  -----
210  ignores config.forceKnownWcs
211  """
212  if self.refObjLoaderrefObjLoader is None:
213  raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
214  import lsstDebug
215  debug = lsstDebug.Info(__name__)
216 
217  expMd = self._getExposureMetadata_getExposureMetadata(exposure)
218 
219  sourceSelection = self.sourceSelector.run(sourceCat)
220 
221  self.log.info("Purged %d sources, leaving %d good sources",
222  len(sourceCat) - len(sourceSelection.sourceCat),
223  len(sourceSelection.sourceCat))
224 
225  loadRes = self.refObjLoaderrefObjLoader.loadPixelBox(
226  bbox=expMd.bbox,
227  wcs=expMd.wcs,
228  filterName=expMd.filterName,
229  photoCalib=expMd.photoCalib,
230  epoch=expMd.epoch,
231  )
232 
233  refSelection = self.referenceSelector.run(loadRes.refCat)
234 
235  matchMeta = self.refObjLoaderrefObjLoader.getMetadataBox(
236  bbox=expMd.bbox,
237  wcs=expMd.wcs,
238  filterName=expMd.filterName,
239  photoCalib=expMd.photoCalib,
240  epoch=expMd.epoch,
241  )
242 
243  if debug.display:
244  frame = int(debug.frame)
246  refCat=refSelection.sourceCat,
247  sourceCat=sourceSelection.sourceCat,
248  exposure=exposure,
249  bbox=expMd.bbox,
250  frame=frame,
251  title="Reference catalog",
252  )
253 
254  res = None
255  wcs = expMd.wcs
256  match_tolerance = None
257  for i in range(self.config.maxIter):
258  iterNum = i + 1
259  try:
260  tryRes = self._matchAndFitWcs_matchAndFitWcs(
261  refCat=refSelection.sourceCat,
262  sourceCat=sourceCat,
263  goodSourceCat=sourceSelection.sourceCat,
264  refFluxField=loadRes.fluxField,
265  bbox=expMd.bbox,
266  wcs=wcs,
267  exposure=exposure,
268  match_tolerance=match_tolerance,
269  )
270  except Exception as e:
271  # if we have had a succeessful iteration then use that; otherwise fail
272  if i > 0:
273  self.log.info("Fit WCS iter %d failed; using previous iteration: %s", iterNum, e)
274  iterNum -= 1
275  break
276  else:
277  raise
278 
279  match_tolerance = tryRes.match_tolerance
280  tryMatchDist = self._computeMatchStatsOnSky_computeMatchStatsOnSky(tryRes.matches)
281  self.log.debug(
282  "Match and fit WCS iteration %d: found %d matches with on-sky distance mean "
283  "= %0.3f +- %0.3f arcsec; max match distance = %0.3f arcsec",
284  iterNum, len(tryRes.matches), tryMatchDist.distMean.asArcseconds(),
285  tryMatchDist.distStdDev.asArcseconds(), tryMatchDist.maxMatchDist.asArcseconds())
286 
287  maxMatchDist = tryMatchDist.maxMatchDist
288  res = tryRes
289  wcs = res.wcs
290  if maxMatchDist.asArcseconds() < self.config.minMatchDistanceArcSec:
291  self.log.debug(
292  "Max match distance = %0.3f arcsec < %0.3f = config.minMatchDistanceArcSec; "
293  "that's good enough",
294  maxMatchDist.asArcseconds(), self.config.minMatchDistanceArcSec)
295  break
296  match_tolerance.maxMatchDist = maxMatchDist
297 
298  self.log.info(
299  "Matched and fit WCS in %d iterations; "
300  "found %d matches with on-sky distance mean and scatter = %0.3f +- %0.3f arcsec",
301  iterNum, len(tryRes.matches), tryMatchDist.distMean.asArcseconds(),
302  tryMatchDist.distStdDev.asArcseconds())
303  if tryMatchDist.distMean.asArcseconds() > self.config.maxMeanDistanceArcsec:
304  raise pipeBase.TaskError(
305  "Fatal astrometry failure detected: mean on-sky distance = %0.3f arcsec > %0.3f "
306  "(maxMeanDistanceArcsec)" %
307  (tryMatchDist.distMean.asArcseconds(), self.config.maxMeanDistanceArcsec))
308  for m in res.matches:
309  if self.usedKeyusedKey:
310  m.second.set(self.usedKeyusedKey, True)
311  exposure.setWcs(res.wcs)
312 
313  # Record the scatter in the exposure metadata
314  md = exposure.getMetadata()
315  md['SFM_ASTROM_OFFSET_MEAN'] = tryMatchDist.distMean.asArcseconds()
316  md['SFM_ASTROM_OFFSET_STD'] = tryMatchDist.distStdDev.asArcseconds()
317 
318  return pipeBase.Struct(
319  refCat=refSelection.sourceCat,
320  matches=res.matches,
321  scatterOnSky=res.scatterOnSky,
322  matchMeta=matchMeta,
323  )
324 
325  @timeMethod
326  def _matchAndFitWcs(self, refCat, sourceCat, goodSourceCat, refFluxField, bbox, wcs, match_tolerance,
327  exposure=None):
328  """Match sources to reference objects and fit a WCS.
329 
330  Parameters
331  ----------
332  refCat : `lsst.afw.table.SimpleCatalog`
333  catalog of reference objects
334  sourceCat : `lsst.afw.table.SourceCatalog`
335  catalog of sources detected on the exposure
336  goodSourceCat : `lsst.afw.table.SourceCatalog`
337  catalog of down-selected good sources detected on the exposure
338  refFluxField : 'str'
339  field of refCat to use for flux
340  bbox : `lsst.geom.Box2I`
341  bounding box of exposure
342  wcs : `lsst.afw.geom.SkyWcs`
343  initial guess for WCS of exposure
344  match_tolerance : `lsst.meas.astrom.MatchTolerance`
345  a MatchTolerance object (or None) specifying
346  internal tolerances to the matcher. See the MatchTolerance
347  definition in the respective matcher for the class definition.
348  exposure : `lsst.afw.image.Exposure`
349  exposure whose WCS is to be fit, or None; used only for the debug
350  display.
351 
352  Returns
353  -------
354  result : `lsst.pipe.base.Struct`
355  Result struct with components:
356 
357  - ``matches``: astrometric matches
358  (`list` of `lsst.afw.table.ReferenceMatch`).
359  - ``wcs``: the fit WCS (lsst.afw.geom.SkyWcs).
360  - ``scatterOnSky`` : median on-sky separation between reference
361  objects and sources in "matches" (`lsst.afw.geom.Angle`).
362  """
363  import lsstDebug
364  debug = lsstDebug.Info(__name__)
365 
366  sourceFluxField = "slot_%sFlux_instFlux" % (self.config.sourceFluxType)
367 
368  matchRes = self.matcher.matchObjectsToSources(
369  refCat=refCat,
370  sourceCat=goodSourceCat,
371  wcs=wcs,
372  sourceFluxField=sourceFluxField,
373  refFluxField=refFluxField,
374  match_tolerance=match_tolerance,
375  )
376  self.log.debug("Found %s matches", len(matchRes.matches))
377  if debug.display:
378  frame = int(debug.frame)
380  refCat=refCat,
381  sourceCat=matchRes.usableSourceCat,
382  matches=matchRes.matches,
383  exposure=exposure,
384  bbox=bbox,
385  frame=frame + 1,
386  title="Initial WCS",
387  )
388 
389  if self.config.doMagnitudeOutlierRejection:
390  matches = self._removeMagnitudeOutliers_removeMagnitudeOutliers(sourceFluxField, refFluxField, matchRes.matches)
391  else:
392  matches = matchRes.matches
393 
394  self.log.debug("Fitting WCS")
395  fitRes = self.wcsFitter.fitWcs(
396  matches=matches,
397  initWcs=wcs,
398  bbox=bbox,
399  refCat=refCat,
400  sourceCat=sourceCat,
401  exposure=exposure,
402  )
403  fitWcs = fitRes.wcs
404  scatterOnSky = fitRes.scatterOnSky
405  if debug.display:
406  frame = int(debug.frame)
408  refCat=refCat,
409  sourceCat=matchRes.usableSourceCat,
410  matches=matches,
411  exposure=exposure,
412  bbox=bbox,
413  frame=frame + 2,
414  title="Fit TAN-SIP WCS",
415  )
416 
417  return pipeBase.Struct(
418  matches=matches,
419  wcs=fitWcs,
420  scatterOnSky=scatterOnSky,
421  match_tolerance=matchRes.match_tolerance,
422  )
423 
424  def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn):
425  """Remove magnitude outliers, computing a simple zeropoint.
426 
427  Parameters
428  ----------
429  sourceFluxField : `str`
430  Field in source catalog for instrumental fluxes.
431  refFluxField : `str`
432  Field in reference catalog for fluxes (nJy).
433  matchesIn : `list` [`lsst.afw.table.ReferenceMatch`]
434  List of source/reference matches input
435 
436  Returns
437  -------
438  matchesOut : `list` [`lsst.afw.table.ReferenceMatch`]
439  List of source/reference matches with magnitude
440  outliers removed.
441  """
442  nMatch = len(matchesIn)
443  sourceMag = np.zeros(nMatch)
444  refMag = np.zeros(nMatch)
445  for i, match in enumerate(matchesIn):
446  sourceMag[i] = -2.5*np.log10(match[1][sourceFluxField])
447  refMag[i] = (match[0][refFluxField]*units.nJy).to_value(units.ABmag)
448 
449  deltaMag = refMag - sourceMag
450  # Protect against negative fluxes and nans in the reference catalog.
451  goodDelta, = np.where(np.isfinite(deltaMag))
452  zp = np.median(deltaMag[goodDelta])
453  # Use median absolute deviation (MAD) for zpSigma.
454  # Also require a minimum scatter to prevent floating-point errors from
455  # rejecting objects in zero-noise tests.
456  zpSigma = np.clip(scipy.stats.median_abs_deviation(deltaMag[goodDelta], scale='normal'),
457  1e-3,
458  None)
459 
460  self.log.info("Rough zeropoint from astrometry matches is %.4f +/- %.4f.",
461  zp, zpSigma)
462 
463  goodStars = goodDelta[(np.abs(deltaMag[goodDelta] - zp)
464  <= self.config.magnitudeOutlierRejectionNSigma*zpSigma)]
465 
466  nOutlier = nMatch - goodStars.size
467  self.log.info("Removed %d magnitude outliers out of %d total astrometry matches.",
468  nOutlier, nMatch)
469 
470  matchesOut = []
471  for matchInd in goodStars:
472  matchesOut.append(matchesIn[matchInd])
473 
474  return matchesOut
def solve(self, exposure, sourceCat)
Definition: astrometry.py:183
def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn)
Definition: astrometry.py:424
def _matchAndFitWcs(self, refCat, sourceCat, goodSourceCat, refFluxField, bbox, wcs, match_tolerance, exposure=None)
Definition: astrometry.py:327
def run(self, sourceCat, exposure)
Definition: astrometry.py:133
def __init__(self, refObjLoader, schema=None, **kwargs)
Definition: astrometry.py:121
def _computeMatchStatsOnSky(self, matchList)
Definition: ref_match.py:208
def loadAndMatch(self, exposure, sourceCat)
Definition: ref_match.py:117
def _getExposureMetadata(self, exposure)
Definition: ref_match.py:237
def displayAstrometry(refCat=None, sourceCat=None, distortedCentroidKey=None, bbox=None, exposure=None, matches=None, frame=1, title="", pause=True)
Definition: display.py:34