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
fitTanSipWcs.py
Go to the documentation of this file.
1 from __future__ import absolute_import, division, print_function
2 
3 import numpy
4 import lsst.afw.image as afwImage
5 import lsst.afw.coord as afwCoord
6 import lsst.afw.geom as afwGeom
7 import lsst.afw.table as afwTable
8 import lsst.pex.config as pexConfig
9 import lsst.pipe.base as pipeBase
10 from .setMatchDistance import setMatchDistance
11 from .sip import makeCreateWcsWithSip
12 
13 __all__ = ["FitTanSipWcsTask", "FitTanSipWcsConfig"]
14 
15 class FitTanSipWcsConfig(pexConfig.Config):
16  order = pexConfig.RangeField(
17  doc = "order of SIP polynomial",
18  dtype = int,
19  default = 4,
20  min = 0,
21  )
22  numIter = pexConfig.RangeField(
23  doc = "number of iterations of fitter (which fits X and Y separately, and so benefits from " + \
24  "a few iterations",
25  dtype = int,
26  default = 3,
27  min = 1,
28  )
29  numRejIter = pexConfig.RangeField(
30  doc = "number of rejection iterations",
31  dtype = int,
32  default = 1,
33  min = 0,
34  )
35  rejSigma = pexConfig.RangeField(
36  doc = "Number of standard deviations for clipping level",
37  dtype = float,
38  default = 3.0,
39  min = 0.0,
40  )
41  maxScatterArcsec = pexConfig.RangeField(
42  doc = "maximum median scatter of a WCS fit beyond which the fit fails (arcsec); " +
43  "be generous, as this is only intended to catch catastrophic failures",
44  dtype = float,
45  default = 10,
46  min = 0,
47  )
48 
49 # The following block adds links to this task from the Task Documentation page.
50 ## \addtogroup LSST_task_documentation
51 ## \{
52 ## \page measAstrom_fitTanSipWcsTask
53 ## \ref FitTanSipWcsTask "FitTanSipWcsTask"
54 ## Fit a TAN-SIP WCS given a list of reference object/source matches
55 ## \}
56 
57 class FitTanSipWcsTask(pipeBase.Task):
58  """!Fit a TAN-SIP WCS given a list of reference object/source matches
59 
60  @anchor FitTanSipWcsTask_
61 
62  @section meas_astrom_fitTanSipWcs_Contents Contents
63 
64  - @ref meas_astrom_fitTanSipWcs_Purpose
65  - @ref meas_astrom_fitTanSipWcs_Initialize
66  - @ref meas_astrom_fitTanSipWcs_IO
67  - @ref meas_astrom_fitTanSipWcs_Schema
68  - @ref meas_astrom_fitTanSipWcs_Config
69  - @ref meas_astrom_fitTanSipWcs_Example
70  - @ref meas_astrom_fitTanSipWcs_Debug
71 
72  @section meas_astrom_fitTanSipWcs_Purpose Description
73 
74  Fit a TAN-SIP WCS given a list of reference object/source matches.
75  See CreateWithSip.h for information about the fitting algorithm.
76 
77  @section meas_astrom_fitTanSipWcs_Initialize Task initialisation
78 
79  @copydoc \_\_init\_\_
80 
81  @section meas_astrom_fitTanSipWcs_IO Invoking the Task
82 
83  @copydoc fitWcs
84 
85  @section meas_astrom_fitTanSipWcs_Config Configuration parameters
86 
87  See @ref FitTanSipWcsConfig
88 
89  @section meas_astrom_fitTanSipWcs_Example A complete example of using FitTanSipWcsTask
90 
91  FitTanSipWcsTask is a subtask of AstrometryTask, which is called by PhotoCalTask.
92  See \ref meas_photocal_photocal_Example.
93 
94  @section meas_astrom_fitTanSipWcs_Debug Debug variables
95 
96  FitTanSipWcsTask does not support any debug variables.
97  """
98  ConfigClass = FitTanSipWcsConfig
99  _DefaultName = "fitWcs"
100 
101  @pipeBase.timeMethod
102  def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None):
103  """!Fit a TAN-SIP WCS from a list of reference object/source matches
104 
105  @param[in,out] matches a list of reference object/source matches
106  (an lsst::afw::table::ReferenceMatchVector)
107  The following fields are read:
108  - match.first (reference object) coord
109  - match.second (source) centroid
110  The following fields are written:
111  - match.first (reference object) centroid,
112  - match.second (source) centroid
113  - match.distance (on sky separation, in radians)
114  @param[in] initWcs initial WCS
115  @param[in] bbox the region over which the WCS will be valid (an lsst:afw::geom::Box2I);
116  if None or an empty box then computed from matches
117  @param[in,out] refCat reference object catalog, or None.
118  If provided then all centroids are updated with the new WCS,
119  otherwise only the centroids for ref objects in matches are updated.
120  Required fields are "centroid_x", "centroid_y", "coord_ra", and "coord_dec".
121  @param[in,out] sourceCat source catalog, or None.
122  If provided then coords are updated with the new WCS;
123  otherwise only the coords for sources in matches are updated.
124  Required fields are "slot_Centroid_x", "slot_Centroid_y", and "coord_ra", and "coord_dec".
125 
126  @return an lsst.pipe.base.Struct with the following fields:
127  - wcs the fit WCS as an lsst.afw.image.Wcs
128  - scatterOnSky median on-sky separation between reference objects and sources in "matches",
129  as an lsst.afw.geom.Angle
130  """
131  if bbox is None:
132  bbox = afwGeom.Box2I()
133 
134  import lsstDebug
135  debug = lsstDebug.Info(__name__)
136 
137  wcs = self.initialWcs(matches, initWcs)
138  rejected = numpy.zeros(len(matches), dtype=bool)
139  for rej in range(self.config.numRejIter):
140  sipObject = self._fitWcs([mm for i, mm in enumerate(matches) if not rejected[i]], wcs)
141  wcs = sipObject.getNewWcs()
142  rejected = self.rejectMatches(matches, wcs, rejected)
143  if rejected.sum() == len(rejected):
144  raise RuntimeError("All matches rejected in iteration %d" % (rej + 1,))
145  if debug.plot:
146  print("Plotting fit after rejection iteration %d/%d" % (rej + 1, self.config.numRejIter))
147  self.plotFit(matches, wcs, rejected)
148  # Final fit after rejection
149  sipObject = self._fitWcs([mm for i, mm in enumerate(matches) if not rejected[i]], wcs)
150  wcs = sipObject.getNewWcs()
151  if debug.plot:
152  print("Plotting final fit")
153  self.plotFit(matches, wcs, rejected)
154 
155  if refCat is not None:
156  self.log.info("Updating centroids in refCat")
157  self.updateRefCentroids(wcs, refList=refCat)
158  else:
159  self.log.warning("Updating reference object centroids in match list; refCat is None")
160  self.updateRefCentroids(wcs, refList=[match.first for match in matches])
161 
162  if sourceCat is not None:
163  self.log.info("Updating coords in sourceCat")
164  self.updateSourceCoords(wcs, sourceList=sourceCat)
165  else:
166  self.log.warning("Updating source coords in match list; sourceCat is None")
167  self.updateSourceCoords(wcs, sourceList=[match.second for match in matches])
168 
169  self.log.info("Updating distance in match list")
170  setMatchDistance(matches)
171 
172  scatterOnSky = sipObject.getScatterOnSky()
173 
174  if scatterOnSky.asArcseconds() > self.config.maxScatterArcsec:
175  raise pipeBase.TaskError(
176  "Fit failed: median scatter on sky = %0.3f arcsec > %0.3f config.maxScatterArcsec" %
177  (scatterOnSky.asArcseconds(), self.config.maxScatterArcsec))
178 
179  return pipeBase.Struct(
180  wcs = wcs,
181  scatterOnSky = scatterOnSky,
182  )
183 
184  def initialWcs(self, matches, wcs):
185  """Generate a guess Wcs from the astrometric matches
186 
187  We create a Wcs anchored at the center of the matches, with the scale
188  of the input Wcs. This is necessary because matching returns only
189  matches with no estimated Wcs, and the input Wcs is a wild guess.
190  We're using the best of each: positions from the matches, and scale
191  from the input Wcs.
192  """
193  crpix = afwGeom.Extent2D(0, 0)
194  crval = afwGeom.Extent3D(0, 0, 0)
195  for mm in matches:
196  crpix += afwGeom.Extent2D(mm.second.getCentroid())
197  crval += afwGeom.Extent3D(mm.first.getCoord().toIcrs().getVector())
198  crpix /= len(matches)
199  crval /= len(matches)
200  newWcs = afwImage.Wcs(afwCoord.IcrsCoord(afwGeom.Point3D(crval)).getPosition(),
201  afwGeom.Point2D(crpix), wcs.getCDMatrix())
202  return newWcs
203 
204  def _fitWcs(self, matches, wcs):
205  """Fit a Wcs based on the matches and a guess Wcs"""
206  for i in range(self.config.numIter):
207  sipObject = makeCreateWcsWithSip(matches, wcs, self.config.order)
208  wcs = sipObject.getNewWcs()
209  return sipObject
210 
211  def rejectMatches(self, matches, wcs, rejected):
212  """Flag deviant matches
213 
214  We return a boolean numpy array indicating whether the corresponding
215  match should be rejected. The previous list of rejections is used
216  so we can calculate uncontaminated statistics.
217  """
218  fit = [wcs.skyToPixel(m.first.getCoord()) for m in matches]
219  dx = numpy.array([ff.getX() - mm.second.getCentroid().getX() for ff, mm in zip(fit, matches)])
220  dy = numpy.array([ff.getY() - mm.second.getCentroid().getY() for ff, mm in zip(fit, matches)])
221  good = numpy.logical_not(rejected)
222  return (dx > self.config.rejSigma*dx[good].std()) | (dy > self.config.rejSigma*dy[good].std())
223 
224  @staticmethod
225  def updateRefCentroids(wcs, refList):
226  """Update centroids in a collection of reference objects, given a WCS
227  """
228  if len(refList) < 1:
229  return
230  schema = refList[0].schema
231  coordKey = afwTable.CoordKey(schema["coord"])
232  centroidKey = afwTable.Point2DKey(schema["centroid"])
233  for refObj in refList:
234  refObj.set(centroidKey, wcs.skyToPixel(refObj.get(coordKey)))
235 
236  @staticmethod
237  def updateSourceCoords(wcs, sourceList):
238  """Update coords in a collection of sources, given a WCS
239  """
240  if len(sourceList) < 1:
241  return
242  schema = sourceList[1].schema
243  srcCoordKey = afwTable.CoordKey(schema["coord"])
244  for src in sourceList:
245  src.set(srcCoordKey, wcs.pixelToSky(src.getCentroid()))
246 
247  def plotFit(self, matches, wcs, rejected):
248  """Plot the fit
249 
250  We create four plots, for all combinations of (dx, dy) against
251  (x, y). Good points are black, while rejected points are red.
252  """
253  import matplotlib.pyplot as plt
254 
255  fit = [wcs.skyToPixel(m.first.getCoord()) for m in matches]
256  x1 = numpy.array([ff.getX() for ff in fit])
257  y1 = numpy.array([ff.getY() for ff in fit])
258  x2 = numpy.array([m.second.getCentroid().getX() for m in matches])
259  y2 = numpy.array([m.second.getCentroid().getY() for m in matches])
260 
261  dx = x1 - x2
262  dy = y1 - y2
263 
264  good = numpy.logical_not(rejected)
265 
266  figure = plt.figure()
267  axes = figure.add_subplot(2, 2, 1)
268  axes.plot(x2[good], dx[good], 'ko')
269  axes.plot(x2[rejected], dx[rejected], 'ro')
270  axes.set_xlabel("x")
271  axes.set_ylabel("dx")
272 
273  axes = figure.add_subplot(2, 2, 2)
274  axes.plot(x2[good], dy[good], 'ko')
275  axes.plot(x2[rejected], dy[rejected], 'ro')
276  axes.set_xlabel("x")
277  axes.set_ylabel("dy")
278 
279  axes = figure.add_subplot(2, 2, 3)
280  axes.plot(y2[good], dx[good], 'ko')
281  axes.plot(y2[rejected], dx[rejected], 'ro')
282  axes.set_xlabel("y")
283  axes.set_ylabel("dx")
284 
285  axes = figure.add_subplot(2, 2, 4)
286  axes.plot(y2[good], dy[good], 'ko')
287  axes.plot(y2[rejected], dy[rejected], 'ro')
288  axes.set_xlabel("y")
289  axes.set_ylabel("dy")
290 
291  plt.show()
def fitWcs
Fit a TAN-SIP WCS from a list of reference object/source matches.
STL namespace.
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
A coordinate class intended to represent absolute positions.
Definition: PSF.h:39
An integer coordinate rectangle.
Definition: Box.h:53
CreateWcsWithSip< MatchT > makeCreateWcsWithSip(std::vector< MatchT > const &matches, afw::image::Wcs const &linearWcs, int const order, afw::geom::Box2I const &bbox=afw::geom::Box2I(), int const ngrid=0)
Factory function for CreateWcsWithSip.
Fit a TAN-SIP WCS given a list of reference object/source matches.
Definition: fitTanSipWcs.py:57
A class to handle Icrs coordinates (inherits from Coord)
Definition: Coord.h:157
A FunctorKey used to get or set celestial coordiantes from a pair of Angle keys.
Definition: aggregates.h:119