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
fitSipDistortion.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__ = ["FitSipDistortionTask", "FitSipDistortionConfig"]
23 
24 
25 import lsst.sphgeom
26 import lsst.pipe.base
27 import lsst.geom
28 import lsst.afw.image
29 import lsst.afw.geom
30 import lsst.afw.display
31 from lsst.utils.timer import timeMethod
32 
33 from .scaledPolynomialTransformFitter import ScaledPolynomialTransformFitter, OutlierRejectionControl
34 from .sipTransform import SipForwardTransform, SipReverseTransform, makeWcs
35 from .makeMatchStatistics import makeMatchStatisticsInRadians
36 
37 from .setMatchDistance import setMatchDistance
38 
39 
41  """Config for FitSipDistortionTask"""
43  doc="Order of SIP polynomial",
44  dtype=int,
45  default=4,
46  min=0,
47  )
49  doc="Number of rejection iterations",
50  dtype=int,
51  default=3,
52  min=0,
53  )
55  doc="Number of standard deviations for clipping level",
56  dtype=float,
57  default=3.0,
58  min=0.0,
59  )
61  doc="Minimum number of matches to reject when sigma-clipping",
62  dtype=int,
63  default=0
64  )
66  doc="Maximum number of matches to reject when sigma-clipping",
67  dtype=int,
68  default=1
69  )
70  maxScatterArcsec = lsst.pex.config.RangeField(
71  doc="Maximum median scatter of a WCS fit beyond which the fit fails (arcsec); "
72  "be generous, as this is only intended to catch catastrophic failures",
73  dtype=float,
74  default=10,
75  min=0,
76  )
77  refUncertainty = lsst.pex.config.Field(
78  doc="RMS uncertainty in reference catalog positions, in pixels. Will be added "
79  "in quadrature with measured uncertainties in the fit.",
80  dtype=float,
81  default=0.25,
82  )
84  doc="Number of X grid points used to invert the SIP reverse transform.",
85  dtype=int,
86  default=100,
87  )
89  doc="Number of Y grid points used to invert the SIP reverse transform.",
90  dtype=int,
91  default=100,
92  )
93  gridBorder = lsst.pex.config.Field(
94  doc="When setting the gird region, how much to extend the image "
95  "bounding box (in pixels) before transforming it to intermediate "
96  "world coordinates using the initial WCS.",
97  dtype=float,
98  default=50.0,
99  )
100 
101 
102 class FitSipDistortionTask(lsst.pipe.base.Task):
103  """Fit a TAN-SIP WCS given a list of reference object/source matches.
104  """
105  ConfigClass = FitSipDistortionConfig
106  _DefaultName = "fitWcs"
107 
108  def __init__(self, **kwargs):
109  lsst.pipe.base.Task.__init__(self, **kwargs)
110  self.outlierRejectionCtrloutlierRejectionCtrl = OutlierRejectionControl()
111  self.outlierRejectionCtrloutlierRejectionCtrl.nClipMin = self.config.nClipMin
112  self.outlierRejectionCtrloutlierRejectionCtrl.nClipMax = self.config.nClipMax
113  self.outlierRejectionCtrloutlierRejectionCtrl.nSigma = self.config.rejSigma
114 
115  @timeMethod
116  def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None):
117  """Fit a TAN-SIP WCS from a list of reference object/source matches.
118 
119  Parameters
120  ----------
121  matches : `list` of `lsst.afw.table.ReferenceMatch`
122  A sequence of reference object/source matches.
123  The following fields are read:
124  - match.first (reference object) coord
125  - match.second (source) centroid
126 
127  The following fields are written:
128  - match.first (reference object) centroid
129  - match.second (source) centroid
130  - match.distance (on sky separation, in radians)
131 
132  initWcs : `lsst.afw.geom.SkyWcs`
133  An initial WCS whose CD matrix is used as the final CD matrix.
134  bbox : `lsst.geom.Box2I`
135  The region over which the WCS will be valid (PARENT pixel coordinates);
136  if `None` or an empty box then computed from matches
137  refCat : `lsst.afw.table.SimpleCatalog`
138  Reference object catalog, or `None`.
139  If provided then all centroids are updated with the new WCS,
140  otherwise only the centroids for ref objects in matches are updated.
141  Required fields are "centroid_x", "centroid_y", "coord_ra", and "coord_dec".
142  sourceCat : `lsst.afw.table.SourceCatalog`
143  Source catalog, or `None`.
144  If provided then coords are updated with the new WCS;
145  otherwise only the coords for sources in matches are updated.
146  Required input fields are "slot_Centroid_x", "slot_Centroid_y",
147  "slot_Centroid_xErr", "slot_Centroid_yErr", and optionally
148  "slot_Centroid_x_y_Cov". The "coord_ra" and "coord_dec" fields
149  will be updated but are not used as input.
150  exposure : `lsst.afw.image.Exposure`
151  An Exposure or other displayable image on which matches can be
152  overplotted. Ignored (and may be `None`) if display-based debugging
153  is not enabled via lsstDebug.
154 
155  Returns
156  -------
157  An lsst.pipe.base.Struct with the following fields:
158  - wcs : `lsst.afw.geom.SkyWcs`
159  The best-fit WCS.
160  - scatterOnSky : `lsst.geom.Angle`
161  The median on-sky separation between reference objects and
162  sources in "matches", as an `lsst.geom.Angle`
163  """
164  import lsstDebug
165  display = lsstDebug.Info(__name__).display
166  displayFrame = lsstDebug.Info(__name__).frame
167  displayPause = lsstDebug.Info(__name__).pause
168 
169  if bbox is None:
170  bbox = lsst.geom.Box2D()
171  for match in matches:
172  bbox.include(match.second.getCentroid())
173  bbox = lsst.geom.Box2I(bbox)
174 
175  wcs = self.makeInitialWcsmakeInitialWcs(matches, initWcs)
176  cdMatrix = lsst.geom.LinearTransform(wcs.getCdMatrix())
177 
178  # Fit the "reverse" mapping from intermediate world coordinates to
179  # pixels, rejecting outliers. Fitting in this direction first makes it
180  # easier to handle the case where we have uncertainty on source
181  # positions but not reference positions. That's the case we have
182  # right now for purely bookeeeping reasons, and it may be the case we
183  # have in the future when we us Gaia as the reference catalog.
184  revFitter = ScaledPolynomialTransformFitter.fromMatches(self.config.order, matches, wcs,
185  self.config.refUncertainty)
186  revFitter.fit()
187  for nIter in range(self.config.numRejIter):
188  revFitter.updateModel()
189  intrinsicScatter = revFitter.updateIntrinsicScatter()
190  clippedSigma, nRejected = revFitter.rejectOutliers(self.outlierRejectionCtrloutlierRejectionCtrl)
191  self.log.debug(
192  "Iteration %s: intrinsic scatter is %4.3f pixels, "
193  "rejected %d outliers at %3.2f sigma.",
194  nIter+1, intrinsicScatter, nRejected, clippedSigma
195  )
196  if display:
197  displayFrame = self.displaydisplay(revFitter, exposure=exposure, bbox=bbox,
198  frame=displayFrame, displayPause=displayPause)
199  revFitter.fit()
200  revScaledPoly = revFitter.getTransform()
201  # Convert the generic ScaledPolynomialTransform result to SIP form
202  # with given CRPIX and CD (this is an exact conversion, up to
203  # floating-point round-off error)
204  sipReverse = SipReverseTransform.convert(revScaledPoly, wcs.getPixelOrigin(), cdMatrix)
205 
206  # Fit the forward mapping to a grid of points created from the reverse
207  # transform. Because that grid needs to be defined in intermediate
208  # world coordinates, and we don't have a good way to get from pixels to
209  # intermediate world coordinates yet (that's what we're fitting), we'll
210  # first grow the box to make it conservatively large...
211  gridBBoxPix = lsst.geom.Box2D(bbox)
212  gridBBoxPix.grow(self.config.gridBorder)
213  # ...and then we'll transform using just the CRPIX offset and CD matrix
214  # linear transform, which is the TAN-only (no SIP distortion, and
215  # hence approximate) mapping from pixels to intermediate world
216  # coordinates.
217  gridBBoxIwc = lsst.geom.Box2D()
218  for point in gridBBoxPix.getCorners():
219  point -= lsst.geom.Extent2D(wcs.getPixelOrigin())
220  gridBBoxIwc.include(cdMatrix(point))
221  fwdFitter = ScaledPolynomialTransformFitter.fromGrid(self.config.order, gridBBoxIwc,
222  self.config.nGridX, self.config.nGridY,
223  revScaledPoly)
224  fwdFitter.fit()
225  # Convert to SIP forward form.
226  fwdScaledPoly = fwdFitter.getTransform()
227  sipForward = SipForwardTransform.convert(fwdScaledPoly, wcs.getPixelOrigin(), cdMatrix)
228 
229  # Make a new WCS from the SIP transform objects and the CRVAL in the
230  # initial WCS.
231  wcs = makeWcs(sipForward, sipReverse, wcs.getSkyOrigin())
232 
233  if refCat is not None:
234  self.log.debug("Updating centroids in refCat")
235  lsst.afw.table.updateRefCentroids(wcs, refList=refCat)
236  else:
237  self.log.warning("Updating reference object centroids in match list; refCat is None")
238  lsst.afw.table.updateRefCentroids(wcs, refList=[match.first for match in matches])
239 
240  if sourceCat is not None:
241  self.log.debug("Updating coords in sourceCat")
242  lsst.afw.table.updateSourceCoords(wcs, sourceList=sourceCat)
243  else:
244  self.log.warning("Updating source coords in match list; sourceCat is None")
245  lsst.afw.table.updateSourceCoords(wcs, sourceList=[match.second for match in matches])
246 
247  self.log.debug("Updating distance in match list")
248  setMatchDistance(matches)
249 
250  stats = makeMatchStatisticsInRadians(wcs, matches, lsst.afw.math.MEDIAN)
251  scatterOnSky = stats.getValue()*lsst.geom.radians
252 
253  if scatterOnSky.asArcseconds() > self.config.maxScatterArcsec:
254  raise lsst.pipe.base.TaskError(
255  "Fit failed: median scatter on sky = %0.3f arcsec > %0.3f config.maxScatterArcsec" %
256  (scatterOnSky.asArcseconds(), self.config.maxScatterArcsec))
257 
258  return lsst.pipe.base.Struct(
259  wcs=wcs,
260  scatterOnSky=scatterOnSky,
261  )
262 
263  def display(self, revFitter, exposure=None, bbox=None, frame=0, pause=True):
264  """Display positions and outlier status overlaid on an image.
265 
266  This method is called by fitWcs when display debugging is enabled. It
267  always drops into pdb before returning to allow interactive inspection,
268  and hence it should never be called in non-interactive contexts.
269 
270  Parameters
271  ----------
272  revFitter : :cpp:class:`lsst::meas::astrom::ScaledPolynomialTransformFitter`
273  Fitter object initialized with `fromMatches` for fitting a "reverse"
274  distortion: the mapping from intermediate world coordinates to
275  pixels.
276  exposure : :cpp:class:`lsst::afw::image::Exposure`
277  An Exposure or other displayable image on which matches can be
278  overplotted.
279  bbox : :cpp:class:`lsst::afw::geom::Box2I`
280  Bounding box of the region on which matches should be plotted.
281  """
282  data = revFitter.getData()
283  disp = lsst.afw.display.getDisplay(frame=frame)
284  if exposure is not None:
285  disp.mtv(exposure)
286  elif bbox is not None:
287  disp.mtv(exposure=lsst.afw.image.ExposureF(bbox))
288  else:
289  raise TypeError("At least one of 'exposure' and 'bbox' must be provided.")
290  data = revFitter.getData()
291  srcKey = lsst.afw.table.Point2DKey(data.schema["src"])
292  srcErrKey = lsst.afw.table.CovarianceMatrix2fKey(data.schema["src"], ["x", "y"])
293  refKey = lsst.afw.table.Point2DKey(data.schema["initial"])
294  modelKey = lsst.afw.table.Point2DKey(data.schema["model"])
295  rejectedKey = data.schema.find("rejected").key
296  with disp.Buffering():
297  for record in data:
298  colors = ((lsst.afw.display.RED, lsst.afw.display.GREEN)
299  if not record.get(rejectedKey) else
300  (lsst.afw.display.MAGENTA, lsst.afw.display.CYAN))
301  rx, ry = record.get(refKey)
302  disp.dot("x", rx, ry, size=10, ctype=colors[0])
303  mx, my = record.get(modelKey)
304  disp.dot("o", mx, my, size=10, ctype=colors[0])
305  disp.line([(rx, ry), (mx, my)], ctype=colors[0])
306  sx, sy = record.get(srcKey)
307  sErr = record.get(srcErrKey)
308  sEllipse = lsst.afw.geom.Quadrupole(sErr[0, 0], sErr[1, 1], sErr[0, 1])
309  disp.dot(sEllipse, sx, sy, ctype=colors[1])
310  if pause or pause is None: # default is to pause
311  print("Dropping into debugger to allow inspection of display. Type 'continue' when done.")
312  import pdb
313  pdb.set_trace()
314  return frame
315  else:
316  return frame + 1 # increment and return the frame for the next iteration.
317 
318  def makeInitialWcs(self, matches, wcs):
319  """Generate a guess Wcs from the astrometric matches
320 
321  We create a Wcs anchored at the center of the matches, with the scale
322  of the input Wcs. This is necessary because the Wcs may have a very
323  approximation position (as is common with telescoped-generated Wcs).
324  We're using the best of each: positions from the matches, and scale
325  from the input Wcs.
326 
327  Parameters
328  ----------
329  matches : list of :cpp:class:`lsst::afw::table::ReferenceMatch`
330  A sequence of reference object/source matches.
331  The following fields are read:
332 
333  - match.first (reference object) coord
334  - match.second (source) centroid
335 
336  wcs : :cpp:class:`lsst::afw::geom::SkyWcs`
337  An initial WCS whose CD matrix is used as the CD matrix of the
338  result.
339 
340  Returns
341  -------
342  newWcs : `lsst.afw.geom.SkyWcs`
343  A new WCS guess.
344  """
345  crpix = lsst.geom.Extent2D(0, 0)
346  crval = lsst.sphgeom.Vector3d(0, 0, 0)
347  for mm in matches:
348  crpix += lsst.geom.Extent2D(mm.second.getCentroid())
349  crval += mm.first.getCoord().getVector()
350  crpix /= len(matches)
351  crval /= len(matches)
352  cd = wcs.getCdMatrix()
353  newWcs = lsst.afw.geom.makeSkyWcs(crpix=lsst.geom.Point2D(crpix),
354  crval=lsst.geom.SpherePoint(crval),
355  cdMatrix=cd)
356  return newWcs
An ellipse core with quadrupole moments as parameters.
Definition: Quadrupole.h:47
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
An integer coordinate rectangle.
Definition: Box.h:55
A 2D linear coordinate transformation.
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None)
def display(self, revFitter, exposure=None, bbox=None, frame=0, pause=True)
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
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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
afw::math::Statistics makeMatchStatisticsInRadians(afw::geom::SkyWcs const &wcs, std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
Compute statistics of on-sky radial separation for a match list, in radians.
std::shared_ptr< afw::geom::SkyWcs > makeWcs(SipForwardTransform const &sipForward, SipReverseTransform const &sipReverse, geom::SpherePoint const &skyOrigin)
Create a new TAN SIP Wcs from a pair of SIP transforms and the sky origin.