LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
anetBasicAstrometry.py
Go to the documentation of this file.
1 from __future__ import absolute_import, division
2 
3 import math
4 import sys
5 
6 import numpy as np # for isfinite()
7 
8 import lsst.daf.base as dafBase
9 import lsst.pex.logging as pexLog
10 from lsst.pex.config import Field, RangeField, ListField
11 import lsst.pex.exceptions as pexExceptions
12 import lsst.pipe.base as pipeBase
13 import lsst.afw.coord as afwCoord
14 import lsst.afw.geom as afwGeom
15 import lsst.afw.table as afwTable
16 import lsst.afw.image as afwImage
17 import lsst.meas.algorithms.utils as maUtils
18 from .loadAstrometryNetObjects import LoadAstrometryNetObjectsTask, LoadMultiIndexes
19 from .display import displayAstrometry
20 
21 from . import sip as astromSip
22 
23 __all__ = ["InitialAstrometry", "ANetBasicAstrometryConfig", "ANetBasicAstrometryTask"]
24 
25 class InitialAstrometry(object):
26  """
27  Object returned by Astrometry.determineWcs
28 
29  getWcs(): sipWcs or tanWcs
30  getMatches(): sipMatches or tanMatches
31 
32  Other fields are:
33  solveQa (PropertyList)
34  tanWcs (Wcs)
35  tanMatches (MatchList)
36  sipWcs (Wcs)
37  sipMatches (MatchList)
38  refCat (lsst::afw::table::SimpleCatalog)
39  matchMeta (PropertyList)
40  """
41  def __init__(self):
42  self.tanWcs = None
43  self.tanMatches = None
44  self.sipWcs = None
45  self.sipMatches = None
46  self.refCat = None
48  self.solveQa = None
49 
50  def getMatches(self):
51  """
52  Get "sipMatches" -- MatchList using the SIP WCS solution, if it
53  exists, or "tanMatches" -- MatchList using the TAN WCS solution
54  otherwise.
55  """
56  return self.sipMatches or self.tanMatches
57 
58  def getWcs(self):
59  """
60  Returns the SIP WCS, if one was found, or a TAN WCS
61  """
62  return self.sipWcs or self.tanWcs
63 
64  matches = property(getMatches)
65  wcs = property(getWcs)
66 
67  ### "Not very pythonic!" complains Paul.
68  # Consider these methods deprecated; if you want these elements, just
69  # .grab them.
70  def getSipWcs(self):
71  return self.sipWcs
72  def getTanWcs(self):
73  return self.tanWcs
74  def getSipMatches(self):
75  return self.sipMatches
76  def getTanMatches(self):
77  return self.tanMatches
78  def getMatchMetadata(self):
79  return self.matchMeta
80  def getSolveQaMetadata(self):
81  return self.solveQa
82 
83 
84 class ANetBasicAstrometryConfig(LoadAstrometryNetObjectsTask.ConfigClass):
85 
86  maxCpuTime = RangeField(
87  doc = "Maximum CPU time to spend solving, in seconds",
88  dtype = float,
89  default = 0.,
90  min = 0.,
91  )
92  matchThreshold = RangeField(
93  doc = "Matching threshold for Astrometry.net solver (log-odds)",
94  dtype = float,
95  default=math.log(1e12),
96  min=math.log(1e6),
97  )
98  maxStars = RangeField(
99  doc = "Maximum number of stars to use in Astrometry.net solving",
100  dtype = int,
101  default=50,
102  min=10,
103  )
104  useWcsPixelScale = Field(
105  doc = "Use the pixel scale from the input exposure\'s WCS headers?",
106  dtype = bool,
107  default = True,
108  )
109  useWcsRaDecCenter = Field(
110  doc="Use the RA,Dec center information from the input exposure\'s WCS headers?",
111  dtype=bool,
112  default=True,
113  )
114  useWcsParity = Field(
115  doc="Use the parity (flip / handedness) of the image from the input exposure\'s WCS headers?",
116  dtype=bool,
117  default=True,
118  )
119  raDecSearchRadius = RangeField(
120  doc = "When useWcsRaDecCenter=True, this is the radius, in degrees, around the RA,Dec center " +
121  "specified in the input exposure\'s WCS to search for a solution.",
122  dtype = float,
123  default=1.0,
124  min=0.0,
125  )
126  pixelScaleUncertainty = RangeField(
127  doc = "Range of pixel scales, around the value in the WCS header, to search. If the value of this field " +
128  "is X and the nominal scale is S, the range searched will be S/X to S*X",
129  dtype = float,
130  default = 1.1,
131  min=1.001,
132  )
133  catalogMatchDist = RangeField(
134  doc = "Matching radius (arcsec) for matching sources to reference objects",
135  dtype = float,
136  default=1.0,
137  min=0.0,
138  )
139  cleaningParameter = RangeField(
140  doc = "Sigma-clipping parameter in sip/cleanBadPoints.py",
141  dtype = float,
142  default=3.0,
143  min=0.0,
144  )
145  calculateSip = Field(
146  doc = "Compute polynomial SIP distortion terms?",
147  dtype = bool,
148  default=True,
149  )
150  sipOrder = RangeField(
151  doc = "Polynomial order of SIP distortion terms",
152  dtype = int,
153  default=4,
154  min=2,
155  )
156  badFlags = ListField(
157  doc = "List of flags which cause a source to be rejected as bad",
158  dtype = str,
159  default = [
160  "slot_Centroid_flag", # bad centroids
161  "base_PixelFlags_flag_edge",
162  "base_PixelFlags_flag_saturated",
163  "base_PixelFlags_flag_crCenter", # cosmic rays
164  ],
165  )
166  allFluxes = Field(
167  doc = "Retrieve all available fluxes (and errors) from catalog?",
168  dtype = bool,
169  default = True,
170  )
171 
172 
173 class ANetBasicAstrometryTask(pipeBase.Task):
174  """!Basic implemeentation of the astrometry.net astrometrical fitter
175 
176  A higher-level class ANetAstrometryTask takes care of dealing with the fact
177  that the initial WCS is probably only a pure TAN SIP, yet we may have
178  significant distortion and a good estimate for that distortion.
179 
180 
181  About Astrometry.net index files (astrometry_net_data):
182 
183  There are three components of an index file: a list of stars
184  (stored as a star kd-tree), a list of quadrangles of stars ("quad
185  file") and a list of the shapes ("codes") of those quadrangles,
186  stored as a code kd-tree.
187 
188  Each index covers a region of the sky, defined by healpix nside
189  and number, and a range of angular scales. In LSST, we share the
190  list of stars in a part of the sky between multiple indexes. That
191  is, the star kd-tree is shared between multiple indices (quads and
192  code kd-trees). In the astrometry.net code, this is called a
193  "multiindex".
194 
195  It is possible to "unload" and "reload" multiindex (and index)
196  objects. When "unloaded", they consume no FILE or mmap resources.
197 
198  The multiindex object holds the star kd-tree and gives each index
199  object it holds a pointer to it, so it is necessary to
200  multiindex_reload_starkd() before reloading the indices it holds.
201  The multiindex_unload() method, on the other hand, unloads its
202  starkd and unloads each index it holds.
203  """
204  ConfigClass = ANetBasicAstrometryConfig
205  _DefaultName = "aNetBasicAstrometry"
206  def __init__(self,
207  config,
208  andConfig=None,
209  **kwargs):
210  """!Construct an ANetBasicAstrometryTask
211 
212  @param[in] config configuration (an instance of self.ConfigClass)
213  @param[in] andConfig astrometry.net data config (an instance of AstromNetDataConfig, or None);
214  if None then use andConfig.py in the astrometry_net_data product (which must be setup)
215 
216  @throw RuntimeError if andConfig is None and the configuration cannot be found,
217  either because astrometry_net_data is not setup in eups
218  or because the setup version does not include the file "andConfig.py"
219  """
220  pipeBase.Task.__init__(self, config=config, **kwargs)
221  self.config = config
222  # this is not a subtask because it cannot safely be retargeted
223  self.refObjLoader = LoadAstrometryNetObjectsTask(config=self.config, andConfig=andConfig, log=self.log,
224  name="loadAN")
225  self.refObjLoader._readIndexFiles()
226 
227  def memusage(self, prefix=''):
228  # Not logging at DEBUG: do nothing
229  if self.log.getThreshold() > pexLog.Log.DEBUG:
230  return
231  from astrometry.util.ttime import get_memusage
232  mu = get_memusage()
233  ss = []
234  for key in ['VmPeak', 'VmSize', 'VmRSS', 'VmData']:
235  if key in mu:
236  ss.append(key + ': ' + ' '.join(mu[key]))
237  if 'mmaps' in mu:
238  ss.append('Mmaps: %i' % len(mu['mmaps']))
239  if 'mmaps_total' in mu:
240  ss.append('Mmaps: %i kB' % (mu['mmaps_total'] / 1024))
241  self.log.logdebug(prefix + 'Memory: ' + ', '.join(ss))
242 
243  def _getImageParams(self, exposure=None, bbox=None, wcs=None, filterName=None, wcsRequired=True):
244  """Get image parameters
245 
246  @param[in] exposure exposure (an afwImage.Exposure) or None
247  @param[in] bbox bounding box (an afwGeom.Box2I) or None; if None then bbox must be specified
248  @param[in] wcs WCS (an afwImage.Wcs) or None; if None then exposure must be specified
249  @param[in] filterName filter name, a string, or None; if None exposure must be specified
250  @param[in] wcsRequired if True then either wcs must be specified or exposure must contain a wcs;
251  if False then the returned wcs may be None
252  @return these items:
253  - bbox bounding box; guaranteed to be set
254  - wcs WCS if known, else None
255  - filterName filter name if known, else None
256  @throw RuntimeError if bbox cannot be determined, or wcs cannot be determined and wcsRequired True
257  """
258  if exposure is not None:
259  if bbox is None:
260  bbox = exposure.getBBox()
261  self.log.logdebug("Setting bbox = %s from exposure metadata" % (bbox,))
262  if wcs is None:
263  self.log.logdebug("Setting wcs from exposure metadata")
264  wcs = exposure.getWcs()
265  if filterName is None:
266  filterName = exposure.getFilter().getName()
267  self.log.logdebug("Setting filterName = %r from exposure metadata" % (filterName,))
268  if bbox is None:
269  raise RuntimeError("bbox or exposure must be specified")
270  if wcs is None and wcsRequired:
271  raise RuntimeError("wcs or exposure (with a WCS) must be specified")
272  return bbox, wcs, filterName
273 
274  def useKnownWcs(self, sourceCat, wcs=None, exposure=None, filterName=None, bbox=None, calculateSip=None):
275  """!Return an InitialAstrometry object, just like determineWcs,
276  but assuming the given input WCS is correct.
277 
278  This involves searching for reference sources within the WCS
279  area, and matching them to the given 'sourceCat'. If
280  'calculateSip' is set, we will try to compute a TAN-SIP
281  distortion correction.
282 
283  @param[in] sourceCat list of detected sources in this image.
284  @param[in] wcs your known WCS, or None to get from exposure
285  @param[in] exposure the exposure holding metadata for this image;
286  if None then you must specify wcs, filterName and bbox
287  @param[in] filterName string, filter name, eg "i", or None to get from exposure`
288  @param[in] bbox bounding box of image, or None to get from exposure
289  @param[in] calculateSip calculate SIP distortion terms for the WCS? If None
290  then use self.config.calculateSip. To disable WCS fitting set calculateSip=False
291 
292  @note this function is also called by 'determineWcs' (via 'determineWcs2'), since the steps are all
293  the same.
294  """
295  # return value:
296  astrom = InitialAstrometry()
297 
298  if calculateSip is None:
299  calculateSip = self.config.calculateSip
300 
301  bbox, wcs, filterName = self._getImageParams(
302  exposure = exposure,
303  bbox = bbox,
304  wcs = wcs,
305  filterName = filterName,
306  wcsRequired = True,
307  )
308  refCat = self.refObjLoader.loadPixelBox(
309  bbox = bbox,
310  wcs = wcs,
311  filterName = filterName,
312  calib = None,
313  ).refCat
314  astrom.refCat = refCat
315  catids = [src.getId() for src in refCat]
316  uids = set(catids)
317  self.log.logdebug('%i reference sources; %i unique IDs' % (len(catids), len(uids)))
318  matches = self._getMatchList(sourceCat, refCat, wcs)
319  uniq = set([sm.second.getId() for sm in matches])
320  if len(matches) != len(uniq):
321  self.log.warn(('The list of matched stars contains duplicate reference source IDs ' +
322  '(%i sources, %i unique ids)') % (len(matches), len(uniq)))
323  if len(matches) == 0:
324  self.log.warn('No matches found between input sources and reference catalogue.')
325  return astrom
326 
327  self.log.logdebug('%i reference objects match input sources using input WCS' % (len(matches)))
328  astrom.tanMatches = matches
329  astrom.tanWcs = wcs
330 
331  srcids = [s.getId() for s in sourceCat]
332  for m in matches:
333  assert(m.second.getId() in srcids)
334  assert(m.second in sourceCat)
335 
336  if calculateSip:
337  sipwcs,matches = self._calculateSipTerms(wcs, refCat, sourceCat, matches, bbox=bbox)
338  if sipwcs == wcs:
339  self.log.logdebug('Failed to find a SIP WCS better than the initial one.')
340  else:
341  self.log.logdebug('%i reference objects match input sources using SIP WCS' %
342  (len(matches),))
343  astrom.sipWcs = sipwcs
344  astrom.sipMatches = matches
345 
346  wcs = astrom.getWcs()
347  # _getMatchList() modifies the source list RA,Dec coordinates.
348  # Here, we make them consistent with the WCS we are returning.
349  for src in sourceCat:
350  src.updateCoord(wcs)
351  astrom.matchMeta = _createMetadata(bbox, wcs, filterName)
352  return astrom
353 
354  def determineWcs(self,
355  sourceCat,
356  exposure,
357  **kwargs):
358  """Find a WCS solution for the given 'sourceCat' in the given
359  'exposure', getting other parameters from config.
360 
361  Valid kwargs include:
362 
363  'radecCenter', an afw.coord.Coord giving the RA,Dec position
364  of the center of the field. This is used to limit the
365  search done by Astrometry.net (to make it faster and load
366  fewer index files, thereby using less memory). Defaults to
367  the RA,Dec center from the exposure's WCS; turn that off
368  with the boolean kwarg 'useRaDecCenter' or config option
369  'useWcsRaDecCenter'
370 
371  'useRaDecCenter', a boolean. Don't use the RA,Dec center from
372  the exposure's initial WCS.
373 
374  'searchRadius', in degrees, to search for a solution around
375  the given 'radecCenter'; default from config option
376  'raDecSearchRadius'.
377 
378  'useParity': parity is the 'flip' of the image. Knowing it
379  reduces the search space (hence time) for Astrometry.net.
380  The parity can be computed from the exposure's WCS (the
381  sign of the determinant of the CD matrix); this option
382  controls whether we do that or force Astrometry.net to
383  search both parities. Default from config.useWcsParity.
384 
385  'pixelScale': afwGeom.Angle, estimate of the angle-per-pixel
386  (ie, arcseconds per pixel). Defaults to a value derived
387  from the exposure's WCS. If enabled, this value, plus or
388  minus config.pixelScaleUncertainty, will be used to limit
389  Astrometry.net's search.
390 
391  'usePixelScale': boolean. Use the pixel scale to limit
392  Astrometry.net's search? Defaults to config.useWcsPixelScale.
393 
394  'filterName', a string, the filter name of this image. Will
395  be mapped through the 'filterMap' config dictionary to a
396  column name in the astrometry_net_data index FITS files.
397  Defaults to the exposure.getFilter() value.
398 
399  'bbox', bounding box of exposure; defaults to exposure.getBBox()
400 
401  """
402  assert(exposure is not None)
403 
404  margs = kwargs.copy()
405  if not 'searchRadius' in margs:
406  margs.update(searchRadius = self.config.raDecSearchRadius * afwGeom.degrees)
407  if not 'usePixelScale' in margs:
408  margs.update(usePixelScale = self.config.useWcsPixelScale)
409  if not 'useRaDecCenter' in margs:
410  margs.update(useRaDecCenter = self.config.useWcsRaDecCenter)
411  if not 'useParity' in margs:
412  margs.update(useParity = self.config.useWcsParity)
413  margs.update(exposure=exposure)
414  return self.determineWcs2(sourceCat=sourceCat, **margs)
415 
416  def determineWcs2(self, sourceCat, **kwargs):
417  """Get a blind astrometric solution for the given catalog of sources.
418 
419  We need:
420  -the image size;
421  -the filter
422 
423  And if available, we can use:
424  -an initial Wcs estimate;
425  --> RA,Dec center
426  --> pixel scale
427  --> "parity"
428 
429  (all of which are metadata of Exposure).
430 
431  filterName: string
432  imageSize: (W,H) integer tuple/iterable
433  pixelScale: afwGeom::Angle per pixel.
434  radecCenter: afwCoord::Coord
435  """
436  wcs,qa = self.getBlindWcsSolution(sourceCat, **kwargs)
437  kw = {}
438  # Keys passed to useKnownWcs
439  for key in ['exposure', 'bbox', 'filterName']:
440  if key in kwargs:
441  kw[key] = kwargs[key]
442  astrom = self.useKnownWcs(sourceCat, wcs=wcs, **kw)
443  astrom.solveQa = qa
444  astrom.tanWcs = wcs
445  return astrom
446 
447  def getBlindWcsSolution(self, sourceCat,
448  exposure=None,
449  wcs=None,
450  bbox=None,
451  radecCenter=None,
452  searchRadius=None,
453  pixelScale=None,
454  filterName=None,
455  doTrim=False,
456  usePixelScale=True,
457  useRaDecCenter=True,
458  useParity=True,
459  searchRadiusScale=2.):
460  if not useRaDecCenter and radecCenter is not None:
461  raise RuntimeError('radecCenter is set, but useRaDecCenter is False. Make up your mind!')
462  if not usePixelScale and pixelScale is not None:
463  raise RuntimeError('pixelScale is set, but usePixelScale is False. Make up your mind!')
464 
465  bbox, wcs, filterName = self._getImageParams(
466  exposure = exposure,
467  bbox = bbox,
468  wcs = wcs,
469  filterName = filterName,
470  wcsRequired = False,
471  )
472 
473  bboxD = afwGeom.Box2D(bbox)
474  xc, yc = bboxD.getCenter()
475  parity = None
476 
477  if wcs is not None:
478  if pixelScale is None:
479  if usePixelScale:
480  pixelScale = wcs.pixelScale()
481  self.log.logdebug('Setting pixel scale estimate = %.3f from given WCS estimate' %
482  (pixelScale.asArcseconds()))
483 
484  if radecCenter is None:
485  if useRaDecCenter:
486  radecCenter = wcs.pixelToSky(xc, yc)
487  self.log.logdebug(('Setting RA,Dec center estimate = (%.3f, %.3f) from given WCS '
488  + 'estimate, using pixel center = (%.1f, %.1f)') %
489  (radecCenter.getLongitude().asDegrees(),
490  radecCenter.getLatitude().asDegrees(), xc, yc))
491 
492  if searchRadius is None:
493  if useRaDecCenter:
494  assert(pixelScale is not None)
495  pixRadius = math.hypot(*bboxD.getDimensions()) / 2
496  searchRadius = (pixelScale * pixRadius * searchRadiusScale)
497  self.log.logdebug(('Using RA,Dec search radius = %.3f deg, from pixel scale, '
498  + 'image size, and searchRadiusScale = %g') %
499  (searchRadius, searchRadiusScale))
500  if useParity:
501  parity = wcs.isFlipped()
502  self.log.logdebug('Using parity = %s' % (parity and 'True' or 'False'))
503 
504  if doTrim:
505  n = len(sourceCat)
506  if exposure is not None:
507  exposureBBoxD = afwGeom.Box2D(exposure.getMaskedImage().getBBox())
508  else:
509  exposureBBoxD = bboxD
510  sourceCat = self._trimBadPoints(sourceCat, exposureBBoxD)
511  self.log.logdebug("Trimming: kept %i of %i sources" % (n, len(sourceCat)))
512 
513  wcs,qa = self._solve(
514  sourceCat = sourceCat,
515  wcs = wcs,
516  bbox = bbox,
517  pixelScale = pixelScale,
518  radecCenter = radecCenter,
519  searchRadius = searchRadius,
520  parity = parity,
521  filterName = filterName,
522  )
523  if wcs is None:
524  raise RuntimeError("Unable to match sources with catalog.")
525  self.log.info('Got astrometric solution from Astrometry.net')
526 
527  rdc = wcs.pixelToSky(xc, yc)
528  self.log.logdebug('New WCS says image center pixel (%.1f, %.1f) -> RA,Dec (%.3f, %.3f)' %
529  (xc, yc, rdc.getLongitude().asDegrees(), rdc.getLatitude().asDegrees()))
530  return wcs,qa
531 
532  def getSipWcsFromWcs(self, wcs, bbox, ngrid=20, linearizeAtCenter=True):
533  """!Get a TAN-SIP WCS, starting from an existing WCS.
534 
535  It uses your WCS to compute a fake grid of corresponding "stars" in pixel and sky coords,
536  and feeds that to the regular SIP code.
537 
538  @param[in] wcs initial WCS
539  @param[in] bbox bounding box of image
540  @param[in] ngrid number of grid points along x and y for fitting (fit at ngrid^2 points)
541  @param[in] linearizeCenter if True, get a linear approximation of the input
542  WCS at the image center and use that as the TAN initialization for
543  the TAN-SIP solution. You probably want this if your WCS has its
544  CRPIX outside the image bounding box.
545  """
546  # Ugh, build src and ref tables
547  srcSchema = afwTable.SourceTable.makeMinimalSchema()
548  key = srcSchema.addField("centroid", type="PointD")
549  srcTable = afwTable.SourceTable.make(srcSchema)
550  srcTable.defineCentroid("centroid")
551  srcs = srcTable
552  refs = afwTable.SimpleTable.make(afwTable.SimpleTable.makeMinimalSchema())
553  cref = []
554  csrc = []
555  (W,H) = bbox.getDimensions()
556  x0, y0 = bbox.getMin()
557  for xx in np.linspace(0., W, ngrid):
558  for yy in np.linspace(0, H, ngrid):
559  src = srcs.makeRecord()
560  src.set(key.getX(), x0 + xx)
561  src.set(key.getY(), y0 + yy)
562  csrc.append(src)
563  rd = wcs.pixelToSky(afwGeom.Point2D(xx + x0, yy + y0))
564  ref = refs.makeRecord()
565  ref.setCoord(rd)
566  cref.append(ref)
567 
568  if linearizeAtCenter:
569  # Linearize the original WCS around the image center to create a
570  # TAN WCS.
571  # Reference pixel in LSST coords
572  crpix = afwGeom.Box2D(bbox).getCenter()
573  crval = wcs.pixelToSky(crpix)
574  crval = crval.getPosition(afwGeom.degrees)
575  # Linearize *AT* crval to get effective CD at crval.
576  # (we use the default skyUnit of degrees as per WCS standard)
577  aff = wcs.linearizePixelToSky(crval)
578  cd = aff.getLinear().getMatrix()
579  wcs = afwImage.Wcs(crval, crpix, cd)
580 
581  return self.getSipWcsFromCorrespondences(wcs, cref, csrc, (W,H),
582  x0=x0, y0=y0)
583 
584 
585  def getSipWcsFromCorrespondences(self, origWcs, refCat, sourceCat, bbox):
586  """Produce a SIP solution given a list of known correspondences.
587 
588  Unlike _calculateSipTerms, this does not iterate the solution;
589  it assumes you have given it a good sets of corresponding stars.
590 
591  NOTE that "refCat" and "sourceCat" are assumed to be the same length;
592  entries "refCat[i]" and "sourceCat[i]" are assumed to be correspondences.
593 
594  @param[in] origWcs the WCS to linearize in order to get the TAN part of the TAN-SIP WCS.
595  @param[in] refCat reference source catalog
596  @param[in] sourceCat source catalog
597  @param[in] bbox bounding box of image
598  """
599  sipOrder = self.config.sipOrder
600  matches = []
601  for ci,si in zip(refCat, sourceCat):
602  matches.append(afwTable.ReferenceMatch(ci, si, 0.))
603 
604  sipObject = astromSip.makeCreateWcsWithSip(matches, origWcs, sipOrder, bbox)
605  return sipObject.getNewWcs()
606 
607  def _calculateSipTerms(self, origWcs, refCat, sourceCat, matches, bbox):
608  """!Iteratively calculate SIP distortions and regenerate matches based on improved WCS.
609 
610  @param[in] origWcs original WCS object, probably (but not necessarily) a TAN WCS;
611  this is used to set the baseline when determining whether a SIP
612  solution is any better; it will be returned if no better SIP solution
613  can be found.
614  @param[in] refCat reference source catalog
615  @param[in] sourceCat sources in the image to be solved
616  @param[in] matches list of supposedly matched sources, using the "origWcs".
617  @param[in] bbox bounding box of image, which is used when finding reverse SIP coefficients.
618  """
619  sipOrder = self.config.sipOrder
620  wcs = origWcs
621 
622  i=0
623  lastScatPix = None
624  while True:
625  try:
626  sipObject = astromSip.makeCreateWcsWithSip(matches, wcs, sipOrder, bbox)
627  if lastScatPix is None:
628  lastScatPix = sipObject.getLinearScatterInPixels()
629  proposedWcs = sipObject.getNewWcs()
630  scatPix = sipObject.getScatterInPixels()
631  self.plotSolution(matches, proposedWcs, bbox.getDimensions())
632  except pexExceptions.Exception as e:
633  self.log.warn('Failed to calculate distortion terms. Error: ' + str(e))
634  break
635 
636  matchSize = len(matches)
637  # use new WCS to get new matchlist.
638  proposedMatchlist = self._getMatchList(sourceCat, refCat, proposedWcs)
639 
640  self.log.logdebug('SIP iteration %i: %i objects match. Median scatter is %g arcsec = %g pixels (vs previous: %i matches, %g pixels)' %
641  (i, len(proposedMatchlist), sipObject.getScatterOnSky().asArcseconds(),
642  scatPix, matchSize, lastScatPix))
643  # Hack convergence tests
644  if len(proposedMatchlist) < matchSize:
645  self.log.logdebug(
646  "Fit WCS: use iter %s because it had more matches than the next iter: %s vs. %s" % \
647  (i, matchSize, len(proposedMatchlist)))
648  break
649  if len(proposedMatchlist) == matchSize and scatPix >= lastScatPix:
650  self.log.logdebug(
651  "Fit WCS: use iter %s because it had less linear scatter than the next iter: %g vs. %g pixels" % \
652  (i, lastScatPix, scatPix))
653  break
654 
655  wcs = proposedWcs
656  matches = proposedMatchlist
657  lastScatPix = scatPix
658  matchSize = len(matches)
659  i += 1
660 
661  return wcs, matches
662 
663  def plotSolution(self, matches, wcs, imageSize):
664  """Plot the solution, when debugging is turned on.
665 
666  @param matches The list of matches
667  @param wcs The Wcs
668  @param imageSize 2-tuple with the image size (W,H)
669  """
670  import lsstDebug
671  display = lsstDebug.Info(__name__).display
672  if not display:
673  return
674 
675  try:
676  import matplotlib.pyplot as plt
677  import numpy
678  except ImportError:
679  print >> sys.stderr, "Unable to import matplotlib"
680  return
681 
682  fig = plt.figure(1)
683  fig.clf()
684  try:
685  fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
686  except: # protect against API changes
687  pass
688 
689  num = len(matches)
690  x = numpy.zeros(num)
691  y = numpy.zeros(num)
692  dx = numpy.zeros(num)
693  dy = numpy.zeros(num)
694  for i, m in enumerate(matches):
695  x[i] = m.second.getX()
696  y[i] = m.second.getY()
697  pixel = wcs.skyToPixel(m.first.getCoord())
698  dx[i] = x[i] - pixel.getX()
699  dy[i] = y[i] - pixel.getY()
700 
701  subplots = maUtils.makeSubplots(fig, 2, 2, xgutter=0.1, ygutter=0.1, pygutter=0.04)
702 
703  def plotNext(x, y, xLabel, yLabel, xMax):
704  ax = subplots.next()
705  ax.set_autoscalex_on(False)
706  ax.set_xbound(lower=0, upper=xMax)
707  ax.scatter(x, y)
708  ax.set_xlabel(xLabel)
709  ax.set_ylabel(yLabel)
710  ax.axhline(0.0)
711 
712  plotNext(x, dx, "x", "dx", imageSize[0])
713  plotNext(x, dy, "x", "dy", imageSize[0])
714  plotNext(y, dx, "y", "dx", imageSize[1])
715  plotNext(y, dy, "y", "dy", imageSize[1])
716 
717  fig.show()
718 
719  while True:
720  try:
721  reply = raw_input("Pausing for inspection, enter to continue... [hpQ] ").strip()
722  except EOFError:
723  reply = "n"
724 
725  reply = reply.split()
726  if reply:
727  reply = reply[0]
728  else:
729  reply = ""
730 
731  if reply in ("", "h", "p", "Q"):
732  if reply == "h":
733  print "h[elp] p[db] Q[uit]"
734  continue
735  elif reply == "p":
736  import pdb; pdb.set_trace()
737  elif reply == "Q":
738  sys.exit(1)
739  break
740 
741  def _getMatchList(self, sourceCat, refCat, wcs):
742  dist = self.config.catalogMatchDist * afwGeom.arcseconds
743  clean = self.config.cleaningParameter
744  matcher = astromSip.MatchSrcToCatalogue(refCat, sourceCat, wcs, dist)
745  matches = matcher.getMatches()
746  if matches is None:
747  # Produce debugging stats...
748  X = [src.getX() for src in sourceCat]
749  Y = [src.getY() for src in sourceCat]
750  R1 = [src.getRa().asDegrees() for src in sourceCat]
751  D1 = [src.getDec().asDegrees() for src in sourceCat]
752  R2 = [src.getRa().asDegrees() for src in refCat]
753  D2 = [src.getDec().asDegrees() for src in refCat]
754  # for src in sourceCat:
755  #self.log.logdebug("source: x,y (%.1f, %.1f), RA,Dec (%.3f, %.3f)" %
756  #(src.getX(), src.getY(), src.getRa().asDegrees(), src.getDec().asDegrees()))
757  #for src in refCat:
758  #self.log.logdebug("ref: RA,Dec (%.3f, %.3f)" %
759  #(src.getRa().asDegrees(), src.getDec().asDegrees()))
760  self.loginfo('_getMatchList: %i sources, %i reference sources' % (len(sourceCat), len(refCat)))
761  if len(sourceCat):
762  self.loginfo(
763  'Source range: x [%.1f, %.1f], y [%.1f, %.1f], RA [%.3f, %.3f], Dec [%.3f, %.3f]' %
764  (min(X), max(X), min(Y), max(Y), min(R1), max(R1), min(D1), max(D1)))
765  if len(refCat):
766  self.loginfo('Reference range: RA [%.3f, %.3f], Dec [%.3f, %.3f]' %
767  (min(R2), max(R2), min(D2), max(D2)))
768  raise RuntimeError('No matches found between image and catalogue')
769  matches = astromSip.cleanBadPoints.clean(matches, wcs, nsigma=clean)
770  return matches
771 
772  def getColumnName(self, filterName, columnMap, default=None):
773  """
774  Returns the column name in the astrometry_net_data index file that will be used
775  for the given filter name.
776 
777  @param filterName Name of filter used in exposure
778  @param columnMap Dict that maps filter names to column names
779  @param default Default column name
780  """
781  filterName = self.config.filterMap.get(filterName, filterName) # Exposure filter --> desired filter
782  try:
783  return columnMap[filterName] # Desired filter --> a_n_d column name
784  except KeyError:
785  self.log.warn("No column in configuration for filter '%s'; using default '%s'" %
786  (filterName, default))
787  return default
788 
789  def _solve(self, sourceCat, wcs, bbox, pixelScale, radecCenter,
790  searchRadius, parity, filterName=None):
791  solver = self.refObjLoader._getSolver()
792 
793  imageSize = bbox.getDimensions()
794  x0, y0 = bbox.getMin()
795 
796  # select sources with valid x,y, flux
797  xybb = afwGeom.Box2D()
798  goodsources = afwTable.SourceCatalog(sourceCat.table)
799  badkeys = [goodsources.getSchema().find(name).key for name in self.config.badFlags]
800 
801  for s in sourceCat:
802  if np.isfinite(s.getX()) and np.isfinite(s.getY()) and np.isfinite(s.getPsfFlux()) \
803  and self._isGoodSource(s, badkeys):
804  goodsources.append(s)
805  xybb.include(afwGeom.Point2D(s.getX() - x0, s.getY() - y0))
806  self.log.info("Number of selected sources for astrometry : %d" %(len(goodsources)))
807  if len(goodsources) < len(sourceCat):
808  self.log.logdebug('Keeping %i of %i sources with finite X,Y positions and PSF flux' %
809  (len(goodsources), len(sourceCat)))
810  self.log.logdebug(('Feeding sources in range x=[%.1f, %.1f], y=[%.1f, %.1f] ' +
811  '(after subtracting x0,y0 = %.1f,%.1f) to Astrometry.net') %
812  (xybb.getMinX(), xybb.getMaxX(), xybb.getMinY(), xybb.getMaxY(), x0, y0))
813  # setStars sorts them by PSF flux.
814  solver.setStars(goodsources, x0, y0)
815  solver.setMaxStars(self.config.maxStars)
816  solver.setImageSize(*imageSize)
817  solver.setMatchThreshold(self.config.matchThreshold)
818  raDecRadius = None
819  if radecCenter is not None:
820  raDecRadius = (radecCenter.getLongitude().asDegrees(),
821  radecCenter.getLatitude().asDegrees(),
822  searchRadius.asDegrees())
823  solver.setRaDecRadius(*raDecRadius)
824  self.log.logdebug('Searching for match around RA,Dec = (%g, %g) with radius %g deg' %
825  raDecRadius)
826 
827  if pixelScale is not None:
828  dscale = self.config.pixelScaleUncertainty
829  scale = pixelScale.asArcseconds()
830  lo = scale / dscale
831  hi = scale * dscale
832  solver.setPixelScaleRange(lo, hi)
833  self.log.logdebug(
834  'Searching for matches with pixel scale = %g +- %g %% -> range [%g, %g] arcsec/pix' %
835  (scale, 100.*(dscale-1.), lo, hi))
836 
837  if parity is not None:
838  solver.setParity(parity)
839  self.log.logdebug('Searching for match with parity = ' + str(parity))
840 
841  # Find and load index files within RA,Dec range and scale range.
842  if radecCenter is not None:
843  multiInds = self.refObjLoader._getMIndexesWithinRange(radecCenter, searchRadius)
844  else:
845  multiInds = self.refObjLoader.multiInds
846  qlo,qhi = solver.getQuadSizeLow(), solver.getQuadSizeHigh()
847 
848  toload_multiInds = set()
849  toload_inds = []
850  for mi in multiInds:
851  for i in range(len(mi)):
852  ind = mi[i]
853  if not ind.overlapsScaleRange(qlo, qhi):
854  continue
855  toload_multiInds.add(mi)
856  toload_inds.append(ind)
857 
858  import lsstDebug
859  if lsstDebug.Info(__name__).display:
860  # Use separate context for display, since astrometry.net can segfault if we don't...
861  with LoadMultiIndexes(toload_multiInds):
862  displayAstrometry(refCat=self.refObjLoader.loadPixelBox(bbox, wcs, filterName).refCat,
863  frame=lsstDebug.Info(__name__).frame, pause=lsstDebug.Info(__name__).pause)
864 
865  with LoadMultiIndexes(toload_multiInds):
866  solver.addIndices(toload_inds)
867  self.memusage('Index files loaded: ')
868 
869  cpulimit = self.config.maxCpuTime
870  solver.run(cpulimit)
871 
872  self.memusage('Solving finished: ')
873 
874  self.memusage('Index files unloaded: ')
875 
876  if solver.didSolve():
877  self.log.logdebug('Solved!')
878  wcs = solver.getWcs()
879  self.log.logdebug('WCS: %s' % wcs.getFitsMetadata().toString())
880 
881  if x0 != 0 or y0 != 0:
882  wcs.shiftReferencePixel(x0, y0)
883  self.log.logdebug('After shifting reference pixel by x0,y0 = (%i,%i), WCS is: %s' %
884  (x0, y0, wcs.getFitsMetadata().toString()))
885 
886  else:
887  self.log.warn('Did not get an astrometric solution from Astrometry.net')
888  wcs = None
889  # Gather debugging info...
890 
891  # -are there any reference stars in the proposed search area?
892  # log the number found and discard the results
893  if radecCenter is not None:
894  self.refObjLoader.loadSkyCircle(radecCenter, searchRadius, filterName)
895 
896  qa = solver.getSolveStats()
897  self.log.logdebug('qa: %s' % qa.toString())
898  return wcs, qa
899 
900  def _isGoodSource(self, candsource, keys):
901  for k in keys:
902  if candsource.get(k):
903  return False
904  return True
905 
906  @staticmethod
907  def _trimBadPoints(sourceCat, bbox, wcs=None):
908  """Remove elements from catalog whose xy positions are not within the given bbox.
909 
910  sourceCat: a Catalog of SimpleRecord or SourceRecord objects
911  bbox: an afwImage.Box2D
912  wcs: if not None, will be used to compute the xy positions on-the-fly;
913  this is required when sources actually contains SimpleRecords.
914 
915  Returns:
916  a list of Source objects with xAstrom,yAstrom within the bbox.
917  """
918  keep = type(sourceCat)(sourceCat.table)
919  for s in sourceCat:
920  point = s.getCentroid() if wcs is None else wcs.skyToPixel(s.getCoord())
921  if bbox.contains(point):
922  keep.append(s)
923  return keep
924 
925  def joinMatchListWithCatalog(self, packedMatches, sourceCat):
926  """
927  This function is required to reconstitute a ReferenceMatchVector after being
928  unpersisted. The persisted form of a ReferenceMatchVector is the
929  normalized Catalog of IDs produced by afw.table.packMatches(), with the result of
930  InitialAstrometry.getMatchMetadata() in the associated tables\' metadata.
931 
932  The "live" form of a matchlist has links to
933  the real record objects that are matched; it is "denormalized".
934  This function takes a normalized match catalog, along with the catalog of
935  sources to which the match catalog refers. It fetches the reference
936  sources that are within range, and then denormalizes the matches
937  -- sets the "matches[*].first" and "matches[*].second" entries
938  to point to the sources in the "sourceCat" argument, and to the
939  reference sources fetched from the astrometry_net_data files.
940 
941  @param[in] packedMatches Unpersisted match list (an lsst.afw.table.BaseCatalog).
942  packedMatches.table.getMetadata() must contain the
943  values from InitialAstrometry.getMatchMetadata()
944  @param[in,out] sourceCat Source catalog used for the 'second' side of the matches
945  (an lsst.afw.table.SourceCatalog). As a side effect,
946  the catalog will be sorted by ID.
947 
948  @return An lsst.afw.table.ReferenceMatchVector of denormalized matches.
949  """
950  matchmeta = packedMatches.table.getMetadata()
951  version = matchmeta.getInt('SMATCHV')
952  if version != 1:
953  raise ValueError('SourceMatchVector version number is %i, not 1.' % version)
954  filterName = matchmeta.getString('FILTER').strip()
955  ctrCoord = afwCoord.IcrsCoord(
956  matchmeta.getDouble('RA') * afwGeom.degrees,
957  matchmeta.getDouble('DEC') * afwGeom.degrees,
958  )
959  rad = matchmeta.getDouble('RADIUS') * afwGeom.degrees
960  refCat = self.refObjLoader.loadSkyCircle(ctrCoord, rad, filterName).refCat
961  refCat.sort()
962  sourceCat.sort()
963  return afwTable.unpackMatches(packedMatches, refCat, sourceCat)
964 
965 
966 def _createMetadata(bbox, wcs, filterName):
967  """
968  Create match metadata entries required for regenerating the catalog
969 
970  @param bbox bounding box of image (pixels)
971  @param filterName Name of filter, used for magnitudes
972  @return Metadata
973  """
974  meta = dafBase.PropertyList()
975 
976  bboxD = afwGeom.Box2D(bbox)
977  cx, cy = bboxD.getCenter()
978  radec = wcs.pixelToSky(cx, cy).toIcrs()
979  meta.add('RA', radec.getRa().asDegrees(), 'field center in degrees')
980  meta.add('DEC', radec.getDec().asDegrees(), 'field center in degrees')
981  pixelRadius = math.hypot(*bboxD.getDimensions())/2.0
982  skyRadius = wcs.pixelScale() * pixelRadius
983  meta.add('RADIUS', skyRadius.asDegrees(),
984  'field radius in degrees, approximate')
985  meta.add('SMATCHV', 1, 'SourceMatchVector version number')
986  if filterName is not None:
987  meta.add('FILTER', filterName, 'LSST filter name for tagalong data')
988  return meta
Basic implemeentation of the astrometry.net astrometrical fitter.
Class for storing ordered metadata with comments.
Definition: PropertyList.h:81
def getSipWcs
&quot;Not very pythonic!&quot; complains Paul.
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
def useKnownWcs
Return an InitialAstrometry object, just like determineWcs, but assuming the given input WCS is corre...
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
Lightweight representation of a geometric match between two records.
Definition: fwd.h:90
def _calculateSipTerms
Iteratively calculate SIP distortions and regenerate matches based on improved WCS.
def getSipWcsFromWcs
Get a TAN-SIP WCS, starting from an existing WCS.
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
A class to handle Icrs coordinates (inherits from Coord)
Definition: Coord.h:157
std::vector< Match< typename Cat1::Record, typename Cat2::Record > > unpackMatches(BaseCatalog const &matches, Cat1 const &cat1, Cat2 const &cat2)
Reconstruct a MatchVector from a BaseCatalog representation of the matches and a pair of catalogs...