LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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 .setMatchDistance import setMatchDistance
34 from .sip import makeCreateWcsWithSip
35 
36 
37 class FitTanSipWcsConfig(pexConfig.Config):
38  """Config for FitTanSipWcsTask."""
39  order = pexConfig.RangeField(
40  doc="order of SIP polynomial",
41  dtype=int,
42  default=2,
43  min=0,
44  )
45  numIter = pexConfig.RangeField(
46  doc="number of iterations of fitter (which fits X and Y separately, and so benefits from "
47  "a few iterations",
48  dtype=int,
49  default=3,
50  min=1,
51  )
52  numRejIter = pexConfig.RangeField(
53  doc="number of rejection iterations",
54  dtype=int,
55  default=1,
56  min=0,
57  )
58  rejSigma = pexConfig.RangeField(
59  doc="Number of standard deviations for clipping level",
60  dtype=float,
61  default=3.0,
62  min=0.0,
63  )
64  maxScatterArcsec = pexConfig.RangeField(
65  doc="maximum median scatter of a WCS fit beyond which the fit fails (arcsec); "
66  "be generous, as this is only intended to catch catastrophic failures",
67  dtype=float,
68  default=10,
69  min=0,
70  )
71 
72 
73 class FitTanSipWcsTask(pipeBase.Task):
74  """Fit a TAN-SIP WCS given a list of reference object/source matches.
75  """
76  ConfigClass = FitTanSipWcsConfig
77  _DefaultName = "fitWcs"
78 
79  @pipeBase.timeMethod
80  def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None):
81  """Fit a TAN-SIP WCS from a list of reference object/source matches
82 
83  Parameters
84  ----------
85  matches : `list` of `lsst.afw.table.ReferenceMatch`
86  The following fields are read:
87 
88  - match.first (reference object) coord
89  - match.second (source) centroid
90 
91  The following fields are written:
92 
93  - match.first (reference object) centroid,
94  - match.second (source) centroid
95  - match.distance (on sky separation, in radians)
96 
97  initWcs : `lsst.afw.geom.SkyWcs`
98  initial WCS
99  bbox : `lsst.geom.Box2I`
100  the region over which the WCS will be valid (an lsst:afw::geom::Box2I);
101  if None or an empty box then computed from matches
102  refCat : `lsst.afw.table.SimpleCatalog`
103  reference object catalog, or None.
104  If provided then all centroids are updated with the new WCS,
105  otherwise only the centroids for ref objects in matches are updated.
106  Required fields are "centroid_x", "centroid_y", "coord_ra", and "coord_dec".
107  sourceCat : `lsst.afw.table.SourceCatalog`
108  source catalog, or None.
109  If provided then coords are updated with the new WCS;
110  otherwise only the coords for sources in matches are updated.
111  Required fields are "slot_Centroid_x", "slot_Centroid_y", and "coord_ra", and "coord_dec".
112  exposure : `lsst.afw.image.Exposure`
113  Ignored; present for consistency with FitSipDistortionTask.
114 
115  Returns
116  -------
117  result : `lsst.pipe.base.Struct`
118  with the following fields:
119 
120  - ``wcs`` : the fit WCS (`lsst.afw.geom.SkyWcs`)
121  - ``scatterOnSky`` : median on-sky separation between reference
122  objects and sources in "matches" (`lsst.afw.geom.Angle`)
123  """
124  if bbox is None:
125  bbox = lsst.geom.Box2I()
126 
127  import lsstDebug
128  debug = lsstDebug.Info(__name__)
129 
130  wcs = self.initialWcsinitialWcs(matches, initWcs)
131  rejected = np.zeros(len(matches), dtype=bool)
132  for rej in range(self.config.numRejIter):
133  sipObject = self._fitWcs_fitWcs([mm for i, mm in enumerate(matches) if not rejected[i]], wcs)
134  wcs = sipObject.getNewWcs()
135  rejected = self.rejectMatchesrejectMatches(matches, wcs, rejected)
136  if rejected.sum() == len(rejected):
137  raise RuntimeError("All matches rejected in iteration %d" % (rej + 1,))
138  self.log.debug(
139  "Iteration {0} of astrometry fitting: rejected {1} outliers, "
140  "out of {2} total matches.".format(
141  rej, rejected.sum(), len(rejected)
142  )
143  )
144  if debug.plot:
145  print("Plotting fit after rejection iteration %d/%d" % (rej + 1, self.config.numRejIter))
146  self.plotFitplotFit(matches, wcs, rejected)
147  # Final fit after rejection
148  sipObject = self._fitWcs_fitWcs([mm for i, mm in enumerate(matches) if not rejected[i]], wcs)
149  wcs = sipObject.getNewWcs()
150  if debug.plot:
151  print("Plotting final fit")
152  self.plotFitplotFit(matches, wcs, rejected)
153 
154  if refCat is not None:
155  self.log.debug("Updating centroids in refCat")
156  afwTable.updateRefCentroids(wcs, refList=refCat)
157  else:
158  self.log.warn("Updating reference object centroids in match list; refCat is None")
159  afwTable.updateRefCentroids(wcs, refList=[match.first for match in matches])
160 
161  if sourceCat is not None:
162  self.log.debug("Updating coords in sourceCat")
163  afwTable.updateSourceCoords(wcs, sourceList=sourceCat)
164  else:
165  self.log.warn("Updating source coords in match list; sourceCat is None")
166  afwTable.updateSourceCoords(wcs, sourceList=[match.second for match in matches])
167 
168  self.log.debug("Updating distance in match list")
169  setMatchDistance(matches)
170 
171  scatterOnSky = sipObject.getScatterOnSky()
172 
173  if scatterOnSky.asArcseconds() > self.config.maxScatterArcsec:
174  raise pipeBase.TaskError(
175  "Fit failed: median scatter on sky = %0.3f arcsec > %0.3f config.maxScatterArcsec" %
176  (scatterOnSky.asArcseconds(), self.config.maxScatterArcsec))
177 
178  return pipeBase.Struct(
179  wcs=wcs,
180  scatterOnSky=scatterOnSky,
181  )
182 
183  def initialWcs(self, matches, wcs):
184  """Generate a guess Wcs from the astrometric matches
185 
186  We create a Wcs anchored at the center of the matches, with the scale
187  of the input Wcs. This is necessary because matching returns only
188  matches with no estimated Wcs, and the input Wcs is a wild guess.
189  We're using the best of each: positions from the matches, and scale
190  from the input Wcs.
191 
192  Parameters
193  ----------
194  matches : `list` of `lsst.afw.table.ReferenceMatch`
195  List of sources matched to references.
196  wcs : `lsst.afw.geom.SkyWcs`
197  Current WCS.
198 
199  Returns
200  -------
201  newWcs : `lsst.afw.geom.SkyWcs`
202  Initial WCS guess from estimated crpix and crval.
203  """
204  crpix = lsst.geom.Extent2D(0, 0)
205  crval = lsst.sphgeom.Vector3d(0, 0, 0)
206  for mm in matches:
207  crpix += lsst.geom.Extent2D(mm.second.getCentroid())
208  crval += mm.first.getCoord().getVector()
209  crpix /= len(matches)
210  crval /= len(matches)
211  newWcs = afwGeom.makeSkyWcs(crpix=lsst.geom.Point2D(crpix),
212  crval=lsst.geom.SpherePoint(crval),
213  cdMatrix=wcs.getCdMatrix())
214  return newWcs
215 
216  def _fitWcs(self, matches, wcs):
217  """Fit a Wcs based on the matches and a guess Wcs.
218 
219  Parameters
220  ----------
221  matches : `list` of `lsst.afw.table.ReferenceMatch`
222  List of sources matched to references.
223  wcs : `lsst.afw.geom.SkyWcs`
224  Current WCS.
225 
226  Returns
227  -------
228  sipObject : `lsst.meas.astrom.sip.CreateWcsWithSip`
229  Fitted SIP object.
230  """
231  for i in range(self.config.numIter):
232  sipObject = makeCreateWcsWithSip(matches, wcs, self.config.order)
233  wcs = sipObject.getNewWcs()
234  return sipObject
235 
236  def rejectMatches(self, matches, wcs, rejected):
237  """Flag deviant matches
238 
239  We return a boolean numpy array indicating whether the corresponding
240  match should be rejected. The previous list of rejections is used
241  so we can calculate uncontaminated statistics.
242 
243  Parameters
244  ----------
245  matches : `list` of `lsst.afw.table.ReferenceMatch`
246  List of sources matched to references.
247  wcs : `lsst.afw.geom.SkyWcs`
248  Fitted WCS.
249  rejected : array-like of `bool`
250  Array of matches rejected from the fit. Unused.
251 
252  Returns
253  -------
254  rejectedMatches : `ndarray` of type `bool`
255  Matched objects found to be outside of tolerance.
256  """
257  fit = [wcs.skyToPixel(m.first.getCoord()) for m in matches]
258  dx = np.array([ff.getX() - mm.second.getCentroid().getX() for ff, mm in zip(fit, matches)])
259  dy = np.array([ff.getY() - mm.second.getCentroid().getY() for ff, mm in zip(fit, matches)])
260  good = np.logical_not(rejected)
261  return (dx > self.config.rejSigma*dx[good].std()) | (dy > self.config.rejSigma*dy[good].std())
262 
263  def plotFit(self, matches, wcs, rejected):
264  """Plot the fit
265 
266  We create four plots, for all combinations of (dx, dy) against
267  (x, y). Good points are black, while rejected points are red.
268 
269  Parameters
270  ----------
271  matches : `list` of `lsst.afw.table.ReferenceMatch`
272  List of sources matched to references.
273  wcs : `lsst.afw.geom.SkyWcs`
274  Fitted WCS.
275  rejected : array-like of `bool`
276  Array of matches rejected from the fit.
277  """
278  try:
279  import matplotlib.pyplot as plt
280  except ImportError as e:
281  self.log.warn("Unable to import matplotlib: %s", e)
282  return
283 
284  fit = [wcs.skyToPixel(m.first.getCoord()) for m in matches]
285  x1 = np.array([ff.getX() for ff in fit])
286  y1 = np.array([ff.getY() for ff in fit])
287  x2 = np.array([m.second.getCentroid().getX() for m in matches])
288  y2 = np.array([m.second.getCentroid().getY() for m in matches])
289 
290  dx = x1 - x2
291  dy = y1 - y2
292 
293  good = np.logical_not(rejected)
294 
295  figure = plt.figure()
296  axes = figure.add_subplot(2, 2, 1)
297  axes.plot(x2[good], dx[good], 'ko')
298  axes.plot(x2[rejected], dx[rejected], 'ro')
299  axes.set_xlabel("x")
300  axes.set_ylabel("dx")
301 
302  axes = figure.add_subplot(2, 2, 2)
303  axes.plot(x2[good], dy[good], 'ko')
304  axes.plot(x2[rejected], dy[rejected], 'ro')
305  axes.set_xlabel("x")
306  axes.set_ylabel("dy")
307 
308  axes = figure.add_subplot(2, 2, 3)
309  axes.plot(y2[good], dx[good], 'ko')
310  axes.plot(y2[rejected], dx[rejected], 'ro')
311  axes.set_xlabel("y")
312  axes.set_ylabel("dx")
313 
314  axes = figure.add_subplot(2, 2, 4)
315  axes.plot(y2[good], dy[good], 'ko')
316  axes.plot(y2[rejected], dy[rejected], 'ro')
317  axes.set_xlabel("y")
318  axes.set_ylabel("dy")
319 
320  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:80
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:526
void updateRefCentroids(geom::SkyWcs const &wcs, ReferenceCollection &refList)
Update centroids in a collection of reference objects.
Definition: wcsUtils.cc:73
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList)
Update sky coordinates in a collection of source objects.
Definition: wcsUtils.cc:96
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.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
STL namespace.