LSSTApplications  17.0+103,17.0+11,17.0+61,18.0.0+13,18.0.0+25,18.0.0+5,18.0.0+52,18.0.0-4-g68ffd23,18.1.0-1-g0001055+8,18.1.0-1-g03d53ef+1,18.1.0-1-g1349e88+28,18.1.0-1-g2505f39+22,18.1.0-1-g380d4d4+27,18.1.0-1-g5315e5e+1,18.1.0-1-g5e4b7ea+10,18.1.0-1-g7e8fceb+1,18.1.0-1-g85f8cd4+23,18.1.0-1-g9a6769a+13,18.1.0-1-ga1a4c1a+22,18.1.0-1-gd55f500+17,18.1.0-12-g42eabe8e+10,18.1.0-14-gd04256d+15,18.1.0-16-g430f6a53+1,18.1.0-17-gd2166b6e4,18.1.0-18-gb5d19ff+1,18.1.0-2-gfbf3545+7,18.1.0-2-gfefb8b5+16,18.1.0-3-g52aa583+13,18.1.0-3-g62b5e86+14,18.1.0-3-g8f4a2b1+17,18.1.0-3-g9bc06b8+7,18.1.0-3-gb69f684+9,18.1.0-4-g1ee41a7+1,18.1.0-5-g6dbcb01+13,18.1.0-5-gc286bb7+3,18.1.0-6-g48bdcd3+2,18.1.0-6-gd05e160+9,18.1.0-7-gc4d902b+2,18.1.0-7-gebc0338+8,18.1.0-9-gae7190a+10,w.2019.38
LSSTDataManagementBasePackage
anetAstrometry.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 <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 __all__ = ["ANetAstrometryConfig", "ANetAstrometryTask", "showAstrometry"]
23 
24 import numpy as np
25 
26 import lsstDebug
28 from lsst.afw.cameraGeom import PIXELS, TAN_PIXELS
29 from lsst.afw.table import Point2DKey, CovarianceMatrix2fKey, updateSourceCoords
30 import lsst.geom as geom
31 import lsst.pex.config as pexConfig
32 import lsst.pipe.base as pipeBase
33 from lsst.meas.astrom import displayAstrometry
34 from lsst.meas.astrom.sip import makeCreateWcsWithSip
35 from .anetBasicAstrometry import ANetBasicAstrometryTask
36 
37 
38 class ANetAstrometryConfig(pexConfig.Config):
39  solver = pexConfig.ConfigurableField(
40  target=ANetBasicAstrometryTask,
41  doc="Basic astrometry solver",
42  )
43  forceKnownWcs = pexConfig.Field(dtype=bool, doc=(
44  "Assume that the input image's WCS is correct, without comparing it to any external reality."
45  " (In contrast to using Astrometry.net). NOTE, if you set this, you probably also want to"
46  " un-set 'solver.calculateSip'; otherwise we'll still try to find a TAN-SIP WCS starting "
47  " from the existing WCS"), default=False)
48  rejectThresh = pexConfig.RangeField(dtype=float, default=3.0, doc="Rejection threshold for Wcs fitting",
49  min=0.0, inclusiveMin=False)
50  rejectIter = pexConfig.RangeField(dtype=int, default=3, doc="Rejection iterations for Wcs fitting",
51  min=0)
52 
53  @property
54  def refObjLoader(self):
55  """An alias, for a uniform interface with the standard AstrometryTask"""
56  return self.solver
57 
58  # \addtogroup LSST_task_documentation
59  # \{
60  # \page measAstrom_anetAstrometryTask
61  # \ref ANetAstrometryTask_ "ANetAstrometryTask"
62  # Use astrometry.net to match input sources with a reference catalog and solve for the Wcs
63  # \}
64 
65 
66 class ANetAstrometryTask(pipeBase.Task):
67  r"""!Use astrometry.net to match input sources with a reference catalog and solve for the Wcs
68 
69  @anchor ANetAstrometryTask_
70 
71  The actual matching and solving is done by the 'solver'; this Task
72  serves as a wrapper for taking into account the known optical distortion.
73 
74  \section pipe_tasks_astrometry_Contents Contents
75 
76  - \ref pipe_tasks_astrometry_Purpose
77  - \ref pipe_tasks_astrometry_Initialize
78  - \ref pipe_tasks_astrometry_IO
79  - \ref pipe_tasks_astrometry_Config
80  - \ref pipe_tasks_astrometry_Debug
81  - \ref pipe_tasks_astrometry_Example
82 
83  \section pipe_tasks_astrometry_Purpose Description
84 
85  \copybrief ANetAstrometryTask
86 
87  \section pipe_tasks_astrometry_Initialize Task initialisation
88 
89  \copydoc \_\_init\_\_
90 
91  \section pipe_tasks_astrometry_IO Invoking the Task
92 
93  \copydoc run
94 
95  \section pipe_tasks_astrometry_Config Configuration parameters
96 
97  See \ref ANetAstrometryConfig
98 
99  \section pipe_tasks_astrometry_Debug Debug variables
100 
101  The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a
102  flag \c -d to import \b debug.py from your \c PYTHONPATH;
103  see \ref baseDebug for more about \b debug.py files.
104 
105  The available variables in ANetAstrometryTask are:
106  <DL>
107  <DT> \c display
108  <DD> If True call showAstrometry while iterating ANetAstrometryConfig.rejectIter times,
109  and also after converging; and call displayAstrometry after applying the distortion correction.
110  <DT> \c frame
111  <DD> display frame to use in showAstrometry and displayAstrometry
112  <DT> \c pause
113  <DD> Pause after showAstrometry and displayAstrometry?
114  </DL>
115 
116  \section pipe_tasks_astrometry_Example A complete example of using ANetAstrometryTask
117 
118  See \ref pipe_tasks_photocal_Example.
119 
120  To investigate the \ref pipe_tasks_astrometry_Debug, put something like
121  \code{.py}
122  import lsstDebug
123  def DebugInfo(name):
124  di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
125  if name in ("lsst.pipe.tasks.anetAstrometry", "lsst.pipe.tasks.anetBasicAstrometry"):
126  di.display = 1
127  di.frame = 1
128  di.pause = True
129 
130  return di
131 
132  lsstDebug.Info = DebugInfo
133  \endcode
134  into your debug.py file and run photoCalTask.py with the \c --debug flag.
135  """
136  ConfigClass = ANetAstrometryConfig
137  _DefaultName = "astrometricSolver"
138 
139  def __init__(self, schema, refObjLoader=None, **kwds):
140  r"""!Create the astrometric calibration task. Most arguments are simply passed onto pipe.base.Task.
141 
142  \param schema An lsst::afw::table::Schema used to create the output lsst.afw.table.SourceCatalog
143  \param refObjLoader The AstrometryTask constructor requires a refObjLoader. In order to make this
144  task retargettable for AstrometryTask it needs to take the same arguments. This argument will be
145  ignored since it uses its own internal loader.
146  \param **kwds keyword arguments to be passed to the lsst.pipe.base.task.Task constructor
147 
148  A centroid field "centroid.distorted" (used internally during the Task's operation)
149  will be added to the schema.
150  """
151  pipeBase.Task.__init__(self, **kwds)
152  self.distortedName = "astrom_distorted"
153  self.centroidXKey = schema.addField(self.distortedName + "_x", type="D",
154  doc="centroid distorted for astrometry solver")
155  self.centroidYKey = schema.addField(self.distortedName + "_y", type="D",
156  doc="centroid distorted for astrometry solver")
157  self.centroidXErrKey = schema.addField(self.distortedName + "_xErr", type="F",
158  doc="centroid distorted err for astrometry solver")
159  self.centroidYErrKey = schema.addField(self.distortedName + "_yErr", type="F",
160  doc="centroid distorted err for astrometry solver")
161  self.centroidFlagKey = schema.addField(self.distortedName + "_flag", type="Flag",
162  doc="centroid distorted flag astrometry solver")
164  self.centroidErrKey = CovarianceMatrix2fKey((self.centroidXErrKey, self.centroidYErrKey))
165  # postpone making the solver subtask because it may not be needed and is expensive to create
166  self.solver = None
167 
168  @pipeBase.timeMethod
169  def run(self, exposure, sourceCat):
170  """!Load reference objects, match sources and optionally fit a WCS
171 
172  This is a thin layer around solve or loadAndMatch, depending on config.forceKnownWcs
173 
174  @param[in,out] exposure exposure whose WCS is to be fit
175  The following are read only:
176  - bbox
177  - photoCalib (may be absent)
178  - filter (may be unset)
179  - detector (if wcs is pure tangent; may be absent)
180  The following are updated:
181  - wcs (the initial value is used as an initial guess, and is required)
182  @param[in] sourceCat catalog of sourceCat detected on the exposure (an lsst.afw.table.SourceCatalog)
183  @return an lsst.pipe.base.Struct with these fields:
184  - refCat reference object catalog of objects that overlap the exposure (with some margin)
185  (an lsst::afw::table::SimpleCatalog)
186  - matches astrometric matches, a list of lsst.afw.table.ReferenceMatch
187  - matchMeta metadata about the field (an lsst.daf.base.PropertyList)
188  """
189  if self.config.forceKnownWcs:
190  return self.loadAndMatch(exposure=exposure, sourceCat=sourceCat)
191  else:
192  return self.solve(exposure=exposure, sourceCat=sourceCat)
193 
194  @pipeBase.timeMethod
195  def solve(self, exposure, sourceCat):
196  r"""!Match with reference sources and calculate an astrometric solution
197 
198  \param[in,out] exposure Exposure to calibrate; wcs is updated
199  \param[in] sourceCat catalog of measured sources (an lsst.afw.table.SourceCatalog)
200  \return a pipeBase.Struct with fields:
201  - refCat reference object catalog of objects that overlap the exposure (with some margin)
202  (an lsst::afw::table::SimpleCatalog)
203  - matches astrometric matches, a list of lsst.afw.table.ReferenceMatch
204  - matchMeta metadata about the field (an lsst.daf.base.PropertyList)
205 
206  The reference catalog actually used is up to the implementation
207  of the solver; it will be manifested in the returned matches as
208  a list of lsst.afw.table.ReferenceMatch objects (\em i.e. of lsst.afw.table.Match with
209  \c first being of type lsst.afw.table.SimpleRecord and \c second type lsst.afw.table.SourceRecord ---
210  the reference object and matched object respectively).
211 
212  \note
213  The input sources have the centroid slot moved to a new column "centroid.distorted"
214  which has the positions corrected for any known optical distortion;
215  the 'solver' (which is instantiated in the 'astrometry' member)
216  should therefore simply use the centroids provided by calling
217  afw.table.Source.getCentroid() on the individual source records. This column \em must
218  be present in the sources table.
219 
220  \note ignores config.forceKnownWcs
221  """
222  results = self._astrometry(sourceCat=sourceCat, exposure=exposure)
223 
224  if results.matches:
225  self.refitWcs(sourceCat=sourceCat, exposure=exposure, matches=results.matches)
226 
227  return results
228 
229  @pipeBase.timeMethod
230  def distort(self, sourceCat, exposure):
231  r"""!Calculate distorted source positions
232 
233  CCD images are often affected by optical distortion that makes
234  the astrometric solution higher order than linear. Unfortunately,
235  most (all?) matching algorithms require that the distortion be
236  small or zero, and so it must be removed. We do this by calculating
237  (un-)distorted positions, based on a known optical distortion model
238  in the Ccd.
239 
240  The distortion correction moves sources, so we return the distorted bounding box.
241 
242  \param[in] exposure Exposure to process
243  \param[in,out] sourceCat SourceCatalog; getX() and getY() will be used as inputs,
244  with distorted points in "centroid.distorted" field.
245  \return bounding box of distorted exposure
246  """
247  detector = exposure.getDetector()
248  pixToTanXYTransform = None
249  if detector is None:
250  self.log.warn("No detector associated with exposure; assuming null distortion")
251  else:
252  pixToTanXYTransform = detector.getTransform(PIXELS, TAN_PIXELS)
253 
254  if pixToTanXYTransform is None:
255  self.log.info("Null distortion correction")
256  for s in sourceCat:
257  s.set(self.centroidKey, s.getCentroid())
258  s.set(self.centroidErrKey, s.getCentroidErr())
259  s.set(self.centroidFlagKey, s.getCentroidFlag())
260  return exposure.getBBox()
261 
262  # Distort source positions
263  self.log.info("Applying distortion correction")
264  for s in sourceCat:
265  centroid = pixToTanXYTransform.forwardTransform(s.getCentroid())
266  s.set(self.centroidKey, centroid)
267  s.set(self.centroidErrKey, s.getCentroidErr())
268  s.set(self.centroidFlagKey, s.getCentroidFlag())
269 
270  # Get distorted image size so that astrometry_net does not clip.
271  bboxD = geom.Box2D()
272  for corner in detector.getCorners(TAN_PIXELS):
273  bboxD.include(corner)
274 
275  if lsstDebug.Info(__name__).display:
276  frame = lsstDebug.Info(__name__).frame
277  pause = lsstDebug.Info(__name__).pause
278  displayAstrometry(sourceCat=sourceCat, distortedCentroidKey=self.centroidKey,
279  exposure=exposure, frame=frame, pause=pause)
280 
281  return geom.Box2I(bboxD)
282 
283  @pipeBase.timeMethod
284  def loadAndMatch(self, exposure, sourceCat, bbox=None):
285  """!Load reference objects overlapping an exposure and match to sources detected on that exposure
286 
287  @param[in] exposure exposure whose WCS is to be fit
288  @param[in] sourceCat catalog of sourceCat detected on the exposure (an lsst.afw.table.SourceCatalog)
289  @param[in] bbox bounding box go use for finding reference objects; if None, use exposure's bbox
290 
291  @return an lsst.pipe.base.Struct with these fields:
292  - refCat reference object catalog of objects that overlap the exposure (with some margin)
293  (an lsst::afw::table::SimpleCatalog)
294  - matches astrometric matches, a list of lsst.afw.table.ReferenceMatch
295  - matchMeta metadata about the field (an lsst.daf.base.PropertyList)
296 
297  @note ignores config.forceKnownWcs
298  """
299  bbox = exposure.getBBox()
300  if not self.solver:
301  self.makeSubtask("solver")
302 
303  astrom = self.solver.useKnownWcs(
304  sourceCat=sourceCat,
305  exposure=exposure,
306  bbox=bbox,
307  calculateSip=False,
308  )
309 
310  if astrom is None or astrom.getWcs() is None:
311  raise RuntimeError("Unable to solve astrometry")
312 
313  matches = astrom.getMatches()
314  matchMeta = astrom.getMatchMetadata()
315  if matches is None or len(matches) == 0:
316  raise RuntimeError("No astrometric matches")
317  self.log.info("%d astrometric matches", len(matches))
318 
319  if self._display:
320  frame = lsstDebug.Info(__name__).frame
321  displayAstrometry(exposure=exposure, sourceCat=sourceCat, matches=matches,
322  frame=frame, pause=False)
323 
324  return pipeBase.Struct(
325  refCat=astrom.refCat,
326  matches=matches,
327  matchMeta=matchMeta,
328  )
329 
330  @pipeBase.timeMethod
331  def _astrometry(self, sourceCat, exposure, bbox=None):
332  r"""!Solve astrometry to produce WCS
333 
334  \param[in] sourceCat Sources on exposure, an lsst.afw.table.SourceCatalog
335  \param[in,out] exposure Exposure to process, an lsst.afw.image.ExposureF or D; wcs is updated
336  \param[in] bbox Bounding box, or None to use exposure
337  \return a pipe.base.Struct with fields:
338  - refCat reference object catalog of objects that overlap the exposure (with some margin)
339  (an lsst::afw::table::SimpleCatalog)
340  - matches astrometric matches, a list of lsst.afw.table.ReferenceMatch
341  - matchMeta metadata about the field (an lsst.daf.base.PropertyList)
342  """
343  self.log.info("Solving astrometry")
344  if bbox is None:
345  bbox = exposure.getBBox()
346 
347  if not self.solver:
348  self.makeSubtask("solver")
349 
350  astrom = self.solver.determineWcs(sourceCat=sourceCat, exposure=exposure, bbox=bbox)
351 
352  if astrom is None or astrom.getWcs() is None:
353  raise RuntimeError("Unable to solve astrometry")
354 
355  matches = astrom.getMatches()
356  matchMeta = astrom.getMatchMetadata()
357  if matches is None or len(matches) == 0:
358  raise RuntimeError("No astrometric matches")
359  self.log.info("%d astrometric matches", len(matches))
360 
361  # Note that this is the Wcs for the provided positions, which may be distorted
362  exposure.setWcs(astrom.getWcs())
363 
364  if self._display:
365  frame = lsstDebug.Info(__name__).frame
366  displayAstrometry(exposure=exposure, sourceCat=sourceCat, matches=matches,
367  frame=frame, pause=False)
368 
369  return pipeBase.Struct(
370  refCat=astrom.refCat,
371  matches=matches,
372  matchMeta=matchMeta,
373  )
374 
375  @pipeBase.timeMethod
376  def refitWcs(self, sourceCat, exposure, matches):
377  """!A final Wcs solution after matching and removing distortion
378 
379  Specifically, fitting the non-linear part, since the linear
380  part has been provided by the matching engine.
381 
382  @param sourceCat Sources on exposure, an lsst.afw.table.SourceCatalog
383  @param exposure Exposure of interest, an lsst.afw.image.ExposureF or D
384  @param matches Astrometric matches, as a list of lsst.afw.table.ReferenceMatch
385 
386  @return the resolved-Wcs object, or None if config.solver.calculateSip is False.
387  """
388  sip = None
389  if self.config.solver.calculateSip:
390  self.log.info("Refitting WCS")
391  origMatches = matches
392  wcs = exposure.getWcs()
393 
394  import lsstDebug
395  display = lsstDebug.Info(__name__).display
396  frame = lsstDebug.Info(__name__).frame
397  pause = lsstDebug.Info(__name__).pause
398 
399  def fitWcs(initialWcs, title=None):
400  """!Do the WCS fitting and display of the results"""
401  sip = makeCreateWcsWithSip(matches, initialWcs, self.config.solver.sipOrder)
402  resultWcs = sip.getNewWcs()
403  if display:
404  showAstrometry(exposure, resultWcs, origMatches, matches, frame=frame,
405  title=title, pause=pause)
406  return resultWcs, sip.getScatterOnSky()
407 
408  numRejected = 0
409  try:
410  for i in range(self.config.rejectIter):
411  wcs, scatter = fitWcs(wcs, title="Iteration %d" % i)
412 
413  ref = np.array([wcs.skyToPixel(m.first.getCoord()) for m in matches])
414  src = np.array([m.second.getCentroid() for m in matches])
415  diff = ref - src
416  rms = diff.std()
417  trimmed = []
418  for d, m in zip(diff, matches):
419  if np.all(np.abs(d) < self.config.rejectThresh*rms):
420  trimmed.append(m)
421  else:
422  numRejected += 1
423  if len(matches) == len(trimmed):
424  break
425  matches = trimmed
426 
427  # Final fit after rejection iterations
428  wcs, scatter = fitWcs(wcs, title="Final astrometry")
429 
431  self.log.warn("Unable to fit SIP: %s", e)
432 
433  self.log.info("Astrometric scatter: %f arcsec (%d matches, %d rejected)",
434  scatter.asArcseconds(), len(matches), numRejected)
435  exposure.setWcs(wcs)
436 
437  # Apply WCS to sources
438  updateSourceCoords(wcs, sourceCat)
439  else:
440  self.log.warn("Not calculating a SIP solution; matches may be suspect")
441 
442  if self._display:
443  frame = lsstDebug.Info(__name__).frame
444  displayAstrometry(exposure=exposure, sourceCat=sourceCat, matches=matches,
445  frame=frame, pause=False)
446 
447  return sip
448 
449 
450 def showAstrometry(exposure, wcs, allMatches, useMatches, frame=0, title=None, pause=False):
451  r"""!Show results of astrometry fitting
452 
453  \param exposure Image to display
454  \param wcs Astrometric solution
455  \param allMatches List of all astrometric matches (including rejects)
456  \param useMatches List of used astrometric matches
457  \param frame Frame number for display
458  \param title Title for display
459  \param pause Pause to allow viewing of the display and optional debugging?
460 
461  - Matches are shown in yellow if used in the Wcs solution, otherwise red
462  - +: Detected objects
463  - x: Catalogue objects
464  """
465  import lsst.afw.display as afwDisplay
466  disp = afwDisplay.Display(frame=frame)
467  disp.mtv(exposure, title=title)
468 
469  useIndices = set(m.second.getId() for m in useMatches)
470 
471  radii = []
472  with disp.Buffering():
473  for i, m in enumerate(allMatches):
474  x, y = m.second.getX(), m.second.getY()
475  pix = wcs.skyToPixel(m.first.getCoord())
476 
477  isUsed = m.second.getId() in useIndices
478  if isUsed:
479  radii.append(np.hypot(pix[0] - x, pix[1] - y))
480 
481  color = afwDisplay.YELLOW if isUsed else afwDisplay.RED
482 
483  disp.dot("+", x, y, size=10, ctype=color)
484  disp.dot("x", pix[0], pix[1], size=10, ctype=color)
485 
486  radii = np.array(radii)
487  print("<dr> = %.4g +- %.4g pixels [%d/%d matches]" % (radii.mean(), radii.std(),
488  len(useMatches), len(allMatches)))
489 
490  if pause:
491  import sys
492  while True:
493  try:
494  reply = input("Debugging? [p]db [q]uit; any other key to continue... ").strip()
495  except EOFError:
496  reply = ""
497 
498  if len(reply) > 1:
499  reply = reply[0]
500  if reply == "p":
501  import pdb
502  pdb.set_trace()
503  elif reply == "q":
504  sys.exit(1)
505  else:
506  break
A floating-point coordinate rectangle geometry.
Definition: Box.h:305
Reports attempts to exceed implementation-defined length limits for some classes. ...
Definition: Runtime.h:76
def loadAndMatch(self, exposure, sourceCat, bbox=None)
Load reference objects overlapping an exposure and match to sources detected on that exposure...
def solve(self, exposure, sourceCat)
Match with reference sources and calculate an astrometric solution.
PointKey< double > Point2DKey
Definition: aggregates.h:118
def distort(self, sourceCat, exposure)
Calculate distorted source positions.
def showAstrometry(exposure, wcs, allMatches, useMatches, frame=0, title=None, pause=False)
Show results of astrometry fitting.
daf::base::PropertySet * set
Definition: fits.cc:884
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList)
Update sky coordinates in a collection of source objects.
Definition: wcsUtils.cc:96
def refitWcs(self, sourceCat, exposure, matches)
A final Wcs solution after matching and removing distortion.
def _astrometry(self, sourceCat, exposure, bbox=None)
Solve astrometry to produce WCS.
def run(self, exposure, sourceCat)
Load reference objects, match sources and optionally fit a WCS.
def displayAstrometry(refCat=None, sourceCat=None, distortedCentroidKey=None, bbox=None, exposure=None, matches=None, frame=1, title="", pause=True)
Definition: display.py:34
CreateWcsWithSip< MatchT > makeCreateWcsWithSip(std::vector< MatchT > const &matches, afw::geom::SkyWcs const &linearWcs, int const order, geom::Box2I const &bbox=geom::Box2I(), int const ngrid=0)
Factory function for CreateWcsWithSip.
An integer coordinate rectangle.
Definition: Box.h:54
def __init__(self, schema, refObjLoader=None, kwds)
Create the astrometric calibration task.
Use astrometry.net to match input sources with a reference catalog and solve for the Wcs...
bool strip
Definition: fits.cc:883