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
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 
32 from .scaledPolynomialTransformFitter import ScaledPolynomialTransformFitter, OutlierRejectionControl
33 from .sipTransform import SipForwardTransform, SipReverseTransform, makeWcs
34 from .makeMatchStatistics import makeMatchStatisticsInRadians
35 
36 from .setMatchDistance import setMatchDistance
37 
38 
40  """Config for FitSipDistortionTask"""
42  doc="Order of SIP polynomial",
43  dtype=int,
44  default=4,
45  min=0,
46  )
48  doc="Number of rejection iterations",
49  dtype=int,
50  default=3,
51  min=0,
52  )
54  doc="Number of standard deviations for clipping level",
55  dtype=float,
56  default=3.0,
57  min=0.0,
58  )
60  doc="Minimum number of matches to reject when sigma-clipping",
61  dtype=int,
62  default=0
63  )
65  doc="Maximum number of matches to reject when sigma-clipping",
66  dtype=int,
67  default=1
68  )
69  maxScatterArcsec = lsst.pex.config.RangeField(
70  doc="Maximum median scatter of a WCS fit beyond which the fit fails (arcsec); "
71  "be generous, as this is only intended to catch catastrophic failures",
72  dtype=float,
73  default=10,
74  min=0,
75  )
76  refUncertainty = lsst.pex.config.Field(
77  doc="RMS uncertainty in reference catalog positions, in pixels. Will be added "
78  "in quadrature with measured uncertainties in the fit.",
79  dtype=float,
80  default=0.25,
81  )
83  doc="Number of X grid points used to invert the SIP reverse transform.",
84  dtype=int,
85  default=100,
86  )
88  doc="Number of Y grid points used to invert the SIP reverse transform.",
89  dtype=int,
90  default=100,
91  )
92  gridBorder = lsst.pex.config.Field(
93  doc="When setting the gird region, how much to extend the image "
94  "bounding box (in pixels) before transforming it to intermediate "
95  "world coordinates using the initial WCS.",
96  dtype=float,
97  default=50.0,
98  )
99 
100 
102  """Fit a TAN-SIP WCS given a list of reference object/source matches.
103  """
104  ConfigClass = FitSipDistortionConfig
105  _DefaultName = "fitWcs"
106 
107  def __init__(self, **kwargs):
108  lsst.pipe.base.Task.__init__(self, **kwargs)
109  self.outlierRejectionCtrloutlierRejectionCtrl = OutlierRejectionControl()
110  self.outlierRejectionCtrloutlierRejectionCtrl.nClipMin = self.configconfig.nClipMin
111  self.outlierRejectionCtrloutlierRejectionCtrl.nClipMax = self.configconfig.nClipMax
112  self.outlierRejectionCtrloutlierRejectionCtrl.nSigma = self.configconfig.rejSigma
113 
114  @lsst.pipe.base.timeMethod
115  def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None):
116  """Fit a TAN-SIP WCS from a list of reference object/source matches.
117 
118  Parameters
119  ----------
120  matches : `list` of `lsst.afw.table.ReferenceMatch`
121  A sequence of reference object/source matches.
122  The following fields are read:
123  - match.first (reference object) coord
124  - match.second (source) centroid
125 
126  The following fields are written:
127  - match.first (reference object) centroid
128  - match.second (source) centroid
129  - match.distance (on sky separation, in radians)
130 
131  initWcs : `lsst.afw.geom.SkyWcs`
132  An initial WCS whose CD matrix is used as the final CD matrix.
133  bbox : `lsst.geom.Box2I`
134  The region over which the WCS will be valid (PARENT pixel coordinates);
135  if `None` or an empty box then computed from matches
136  refCat : `lsst.afw.table.SimpleCatalog`
137  Reference object catalog, or `None`.
138  If provided then all centroids are updated with the new WCS,
139  otherwise only the centroids for ref objects in matches are updated.
140  Required fields are "centroid_x", "centroid_y", "coord_ra", and "coord_dec".
141  sourceCat : `lsst.afw.table.SourceCatalog`
142  Source catalog, or `None`.
143  If provided then coords are updated with the new WCS;
144  otherwise only the coords for sources in matches are updated.
145  Required input fields are "slot_Centroid_x", "slot_Centroid_y",
146  "slot_Centroid_xErr", "slot_Centroid_yErr", and optionally
147  "slot_Centroid_x_y_Cov". The "coord_ra" and "coord_dec" fields
148  will be updated but are not used as input.
149  exposure : `lsst.afw.image.Exposure`
150  An Exposure or other displayable image on which matches can be
151  overplotted. Ignored (and may be `None`) if display-based debugging
152  is not enabled via lsstDebug.
153 
154  Returns
155  -------
156  An lsst.pipe.base.Struct with the following fields:
157  - wcs : `lsst.afw.geom.SkyWcs`
158  The best-fit WCS.
159  - scatterOnSky : `lsst.geom.Angle`
160  The median on-sky separation between reference objects and
161  sources in "matches", as an `lsst.geom.Angle`
162  """
163  import lsstDebug
164  display = lsstDebug.Info(__name__).display
165  displayFrame = lsstDebug.Info(__name__).frame
166  displayPause = lsstDebug.Info(__name__).pause
167 
168  if bbox is None:
169  bbox = lsst.geom.Box2D()
170  for match in matches:
171  bbox.include(match.second.getCentroid())
172  bbox = lsst.geom.Box2I(bbox)
173 
174  wcs = self.makeInitialWcsmakeInitialWcs(matches, initWcs)
175  cdMatrix = lsst.geom.LinearTransform(wcs.getCdMatrix())
176 
177  # Fit the "reverse" mapping from intermediate world coordinates to
178  # pixels, rejecting outliers. Fitting in this direction first makes it
179  # easier to handle the case where we have uncertainty on source
180  # positions but not reference positions. That's the case we have
181  # right now for purely bookeeeping reasons, and it may be the case we
182  # have in the future when we us Gaia as the reference catalog.
183  revFitter = ScaledPolynomialTransformFitter.fromMatches(self.configconfig.order, matches, wcs,
184  self.configconfig.refUncertainty)
185  revFitter.fit()
186  for nIter in range(self.configconfig.numRejIter):
187  revFitter.updateModel()
188  intrinsicScatter = revFitter.updateIntrinsicScatter()
189  clippedSigma, nRejected = revFitter.rejectOutliers(self.outlierRejectionCtrloutlierRejectionCtrl)
190  self.loglog.debug(
191  "Iteration {0}: intrinsic scatter is {1:4.3f} pixels, "
192  "rejected {2} outliers at {3:3.2f} sigma.".format(
193  nIter+1, intrinsicScatter, nRejected, clippedSigma
194  )
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.configconfig.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.configconfig.order, gridBBoxIwc,
222  self.configconfig.nGridX, self.configconfig.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.loglog.debug("Updating centroids in refCat")
235  lsst.afw.table.updateRefCentroids(wcs, refList=refCat)
236  else:
237  self.loglog.warn("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.loglog.debug("Updating coords in sourceCat")
242  lsst.afw.table.updateSourceCoords(wcs, sourceList=sourceCat)
243  else:
244  self.loglog.warn("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.loglog.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.configconfig.maxScatterArcsec:
255  "Fit failed: median scatter on sky = %0.3f arcsec > %0.3f config.maxScatterArcsec" %
256  (scatterOnSky.asArcseconds(), self.configconfig.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:526
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:73
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList)
Update sky coordinates in a collection of source objects.
Definition: wcsUtils.cc:96
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.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174