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
fitTanSipWcs.py
Go to the documentation of this file.
1 # This file is part of meas_astrom.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <https://www.gnu.org/licenses/>.
21 
22 __all__ = ["FitTanSipWcsTask", "FitTanSipWcsConfig"]
23 
24 
25 import numpy as np
26 
27 import lsst.geom
28 import lsst.sphgeom
29 import lsst.afw.geom as afwGeom
30 import lsst.afw.table as afwTable
31 import lsst.pex.config as pexConfig
32 import lsst.pipe.base as pipeBase
33 from lsst.utils.timer import timeMethod
34 from .setMatchDistance import setMatchDistance
35 from .sip import makeCreateWcsWithSip
36 
37 
38 class FitTanSipWcsConfig(pexConfig.Config):
39  """Config for FitTanSipWcsTask."""
40  order = pexConfig.RangeField(
41  doc="order of SIP polynomial",
42  dtype=int,
43  default=2,
44  min=0,
45  )
46  numIter = pexConfig.RangeField(
47  doc="number of iterations of fitter (which fits X and Y separately, and so benefits from "
48  "a few iterations",
49  dtype=int,
50  default=3,
51  min=1,
52  )
53  numRejIter = pexConfig.RangeField(
54  doc="number of rejection iterations",
55  dtype=int,
56  default=1,
57  min=0,
58  )
59  rejSigma = pexConfig.RangeField(
60  doc="Number of standard deviations for clipping level",
61  dtype=float,
62  default=3.0,
63  min=0.0,
64  )
65  maxScatterArcsec = pexConfig.RangeField(
66  doc="maximum median scatter of a WCS fit beyond which the fit fails (arcsec); "
67  "be generous, as this is only intended to catch catastrophic failures",
68  dtype=float,
69  default=10,
70  min=0,
71  )
72 
73 
74 class FitTanSipWcsTask(pipeBase.Task):
75  """Fit a TAN-SIP WCS given a list of reference object/source matches.
76  """
77  ConfigClass = FitTanSipWcsConfig
78  _DefaultName = "fitWcs"
79 
80  @timeMethod
81  def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None):
82  """Fit a TAN-SIP WCS from a list of reference object/source matches
83 
84  Parameters
85  ----------
86  matches : `list` of `lsst.afw.table.ReferenceMatch`
87  The following fields are read:
88 
89  - match.first (reference object) coord
90  - match.second (source) centroid
91 
92  The following fields are written:
93 
94  - match.first (reference object) centroid,
95  - match.second (source) centroid
96  - match.distance (on sky separation, in radians)
97 
98  initWcs : `lsst.afw.geom.SkyWcs`
99  initial WCS
100  bbox : `lsst.geom.Box2I`
101  the region over which the WCS will be valid (an lsst:afw::geom::Box2I);
102  if None or an empty box then computed from matches
103  refCat : `lsst.afw.table.SimpleCatalog`
104  reference object catalog, or None.
105  If provided then all centroids are updated with the new WCS,
106  otherwise only the centroids for ref objects in matches are updated.
107  Required fields are "centroid_x", "centroid_y", "coord_ra", and "coord_dec".
108  sourceCat : `lsst.afw.table.SourceCatalog`
109  source catalog, or None.
110  If provided then coords are updated with the new WCS;
111  otherwise only the coords for sources in matches are updated.
112  Required fields are "slot_Centroid_x", "slot_Centroid_y", and "coord_ra", and "coord_dec".
113  exposure : `lsst.afw.image.Exposure`
114  Ignored; present for consistency with FitSipDistortionTask.
115 
116  Returns
117  -------
118  result : `lsst.pipe.base.Struct`
119  with the following fields:
120 
121  - ``wcs`` : the fit WCS (`lsst.afw.geom.SkyWcs`)
122  - ``scatterOnSky`` : median on-sky separation between reference
123  objects and sources in "matches" (`lsst.afw.geom.Angle`)
124  """
125  if bbox is None:
126  bbox = lsst.geom.Box2I()
127 
128  import lsstDebug
129  debug = lsstDebug.Info(__name__)
130 
131  wcs = self.initialWcsinitialWcs(matches, initWcs)
132  rejected = np.zeros(len(matches), dtype=bool)
133  for rej in range(self.config.numRejIter):
134  sipObject = self._fitWcs_fitWcs([mm for i, mm in enumerate(matches) if not rejected[i]], wcs)
135  wcs = sipObject.getNewWcs()
136  rejected = self.rejectMatchesrejectMatches(matches, wcs, rejected)
137  if rejected.sum() == len(rejected):
138  raise RuntimeError("All matches rejected in iteration %d" % (rej + 1,))
139  self.log.debug(
140  "Iteration %d of astrometry fitting: rejected %f outliers, out of %d total matches.",
141  rej, rejected.sum(), len(rejected)
142  )
143  if debug.plot:
144  print("Plotting fit after rejection iteration %d/%d" % (rej + 1, self.config.numRejIter))
145  self.plotFitplotFit(matches, wcs, rejected)
146  # Final fit after rejection
147  sipObject = self._fitWcs_fitWcs([mm for i, mm in enumerate(matches) if not rejected[i]], wcs)
148  wcs = sipObject.getNewWcs()
149  if debug.plot:
150  print("Plotting final fit")
151  self.plotFitplotFit(matches, wcs, rejected)
152 
153  if refCat is not None:
154  self.log.debug("Updating centroids in refCat")
155  afwTable.updateRefCentroids(wcs, refList=refCat)
156  else:
157  self.log.warning("Updating reference object centroids in match list; refCat is None")
158  afwTable.updateRefCentroids(wcs, refList=[match.first for match in matches])
159 
160  if sourceCat is not None:
161  self.log.debug("Updating coords in sourceCat")
162  afwTable.updateSourceCoords(wcs, sourceList=sourceCat)
163  else:
164  self.log.warning("Updating source coords in match list; sourceCat is None")
165  afwTable.updateSourceCoords(wcs, sourceList=[match.second for match in matches])
166 
167  self.log.debug("Updating distance in match list")
168  setMatchDistance(matches)
169 
170  scatterOnSky = sipObject.getScatterOnSky()
171 
172  if scatterOnSky.asArcseconds() > self.config.maxScatterArcsec:
173  raise pipeBase.TaskError(
174  "Fit failed: median scatter on sky = %0.3f arcsec > %0.3f config.maxScatterArcsec" %
175  (scatterOnSky.asArcseconds(), self.config.maxScatterArcsec))
176 
177  return pipeBase.Struct(
178  wcs=wcs,
179  scatterOnSky=scatterOnSky,
180  )
181 
182  def initialWcs(self, matches, wcs):
183  """Generate a guess Wcs from the astrometric matches
184 
185  We create a Wcs anchored at the center of the matches, with the scale
186  of the input Wcs. This is necessary because matching returns only
187  matches with no estimated Wcs, and the input Wcs is a wild guess.
188  We're using the best of each: positions from the matches, and scale
189  from the input Wcs.
190 
191  Parameters
192  ----------
193  matches : `list` of `lsst.afw.table.ReferenceMatch`
194  List of sources matched to references.
195  wcs : `lsst.afw.geom.SkyWcs`
196  Current WCS.
197 
198  Returns
199  -------
200  newWcs : `lsst.afw.geom.SkyWcs`
201  Initial WCS guess from estimated crpix and crval.
202  """
203  crpix = lsst.geom.Extent2D(0, 0)
204  crval = lsst.sphgeom.Vector3d(0, 0, 0)
205  for mm in matches:
206  crpix += lsst.geom.Extent2D(mm.second.getCentroid())
207  crval += mm.first.getCoord().getVector()
208  crpix /= len(matches)
209  crval /= len(matches)
210  newWcs = afwGeom.makeSkyWcs(crpix=lsst.geom.Point2D(crpix),
211  crval=lsst.geom.SpherePoint(crval),
212  cdMatrix=wcs.getCdMatrix())
213  return newWcs
214 
215  def _fitWcs(self, matches, wcs):
216  """Fit a Wcs based on the matches and a guess Wcs.
217 
218  Parameters
219  ----------
220  matches : `list` of `lsst.afw.table.ReferenceMatch`
221  List of sources matched to references.
222  wcs : `lsst.afw.geom.SkyWcs`
223  Current WCS.
224 
225  Returns
226  -------
227  sipObject : `lsst.meas.astrom.sip.CreateWcsWithSip`
228  Fitted SIP object.
229  """
230  for i in range(self.config.numIter):
231  sipObject = makeCreateWcsWithSip(matches, wcs, self.config.order)
232  wcs = sipObject.getNewWcs()
233  return sipObject
234 
235  def rejectMatches(self, matches, wcs, rejected):
236  """Flag deviant matches
237 
238  We return a boolean numpy array indicating whether the corresponding
239  match should be rejected. The previous list of rejections is used
240  so we can calculate uncontaminated statistics.
241 
242  Parameters
243  ----------
244  matches : `list` of `lsst.afw.table.ReferenceMatch`
245  List of sources matched to references.
246  wcs : `lsst.afw.geom.SkyWcs`
247  Fitted WCS.
248  rejected : array-like of `bool`
249  Array of matches rejected from the fit. Unused.
250 
251  Returns
252  -------
253  rejectedMatches : `ndarray` of type `bool`
254  Matched objects found to be outside of tolerance.
255  """
256  fit = [wcs.skyToPixel(m.first.getCoord()) for m in matches]
257  dx = np.array([ff.getX() - mm.second.getCentroid().getX() for ff, mm in zip(fit, matches)])
258  dy = np.array([ff.getY() - mm.second.getCentroid().getY() for ff, mm in zip(fit, matches)])
259  good = np.logical_not(rejected)
260  return (dx > self.config.rejSigma*dx[good].std()) | (dy > self.config.rejSigma*dy[good].std())
261 
262  def plotFit(self, matches, wcs, rejected):
263  """Plot the fit
264 
265  We create four plots, for all combinations of (dx, dy) against
266  (x, y). Good points are black, while rejected points are red.
267 
268  Parameters
269  ----------
270  matches : `list` of `lsst.afw.table.ReferenceMatch`
271  List of sources matched to references.
272  wcs : `lsst.afw.geom.SkyWcs`
273  Fitted WCS.
274  rejected : array-like of `bool`
275  Array of matches rejected from the fit.
276  """
277  try:
278  import matplotlib.pyplot as plt
279  except ImportError as e:
280  self.log.warning("Unable to import matplotlib: %s", e)
281  return
282 
283  fit = [wcs.skyToPixel(m.first.getCoord()) for m in matches]
284  x1 = np.array([ff.getX() for ff in fit])
285  y1 = np.array([ff.getY() for ff in fit])
286  x2 = np.array([m.second.getCentroid().getX() for m in matches])
287  y2 = np.array([m.second.getCentroid().getY() for m in matches])
288 
289  dx = x1 - x2
290  dy = y1 - y2
291 
292  good = np.logical_not(rejected)
293 
294  figure = plt.figure()
295  axes = figure.add_subplot(2, 2, 1)
296  axes.plot(x2[good], dx[good], 'ko')
297  axes.plot(x2[rejected], dx[rejected], 'ro')
298  axes.set_xlabel("x")
299  axes.set_ylabel("dx")
300 
301  axes = figure.add_subplot(2, 2, 2)
302  axes.plot(x2[good], dy[good], 'ko')
303  axes.plot(x2[rejected], dy[rejected], 'ro')
304  axes.set_xlabel("x")
305  axes.set_ylabel("dy")
306 
307  axes = figure.add_subplot(2, 2, 3)
308  axes.plot(y2[good], dx[good], 'ko')
309  axes.plot(y2[rejected], dx[rejected], 'ro')
310  axes.set_xlabel("y")
311  axes.set_ylabel("dx")
312 
313  axes = figure.add_subplot(2, 2, 4)
314  axes.plot(y2[good], dy[good], 'ko')
315  axes.plot(y2[rejected], dy[rejected], 'ro')
316  axes.set_xlabel("y")
317  axes.set_ylabel("dy")
318 
319  plt.show()
An integer coordinate rectangle.
Definition: Box.h:55
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None)
Definition: fitTanSipWcs.py:81
def plotFit(self, matches, wcs, rejected)
def rejectMatches(self, matches, wcs, rejected)
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:521
void updateRefCentroids(geom::SkyWcs const &wcs, ReferenceCollection &refList)
Update centroids in a collection of reference objects.
Definition: wcsUtils.cc:72
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList)
Update sky coordinates in a collection of source objects.
Definition: wcsUtils.cc:95
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.
STL namespace.