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
fitAffineWcs.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__ = ["FitAffineWcsTask", "FitAffineWcsConfig", "TransformedSkyWcsMaker"]
23 
24 
25 import astshim
26 import numpy as np
27 from scipy.optimize import least_squares
28 
29 from lsst.afw.geom import makeSkyWcs, SkyWcs
30 import lsst.afw.math
31 from lsst.geom import Point2D, degrees, arcseconds, radians
32 import lsst.pex.config as pexConfig
33 import lsst.pipe.base as pipeBase
34 from lsst.utils.timer import timeMethod
35 
36 from .makeMatchStatistics import makeMatchStatisticsInRadians
37 from .setMatchDistance import setMatchDistance
38 
39 
40 def _chiFunc(x, refPoints, srcPixels, wcsMaker):
41  """Function to minimize to fit the shift and rotation in the WCS.
42 
43  Parameters
44  ----------
45  x : `numpy.ndarray`
46  Current fit values to test. Float values in array are:
47 
48  - ``bearingTo``: Direction to move the wcs coord in.
49  - ``separation``: Distance along sphere to move wcs coord in.
50  - ``affine0,0``: [0, 0] value of the 2x2 affine transform matrix.
51  - ``affine0,1``: [0, 1] value of the 2x2 affine transform matrix.
52  - ``affine1,0``: [1, 0] value of the 2x2 affine transform matrix.
53  - ``affine1,1``: [1, 1] value of the 2x2 affine transform matrix.
54  refPoints : `list` of `lsst.afw.geom.SpherePoint`
55  Reference object on Sky locations.
56  srcPixels : `list` of `lsst.geom.Point2D`
57  Source object positions on the pixels.
58  wcsMaker : `TransformedSkyWcsMaker`
59  Container class for producing the updated Wcs.
60 
61  Returns
62  -------
63  outputSeparations : `list` of `float`
64  Separation between predicted source location and reference location in
65  radians.
66  """
67  wcs = wcsMaker.makeWcs(x[:2], x[2:].reshape((2, 2)))
68 
69  outputSeparations = []
70  # Fit both sky to pixel and pixel to sky to avoid any non-invertible
71  # affine matrices.
72  for ref, src in zip(refPoints, srcPixels):
73  skySep = ref.getTangentPlaneOffset(wcs.pixelToSky(src))
74  outputSeparations.append(skySep[0].asArcseconds())
75  outputSeparations.append(skySep[1].asArcseconds())
76  xySep = src - wcs.skyToPixel(ref)
77  # Convert the pixel separations to units, arcseconds to match units
78  # of sky separation.
79  outputSeparations.append(
80  xySep[0] * wcs.getPixelScale(src).asArcseconds())
81  outputSeparations.append(
82  xySep[1] * wcs.getPixelScale(src).asArcseconds())
83 
84  return outputSeparations
85 
86 
87 # Keeping this around for now in case any of the fit parameters need to be
88 # configurable. Likely the maximum allowed shift magnitude (parameter 2 in the
89 # fit.)
90 class FitAffineWcsConfig(pexConfig.Config):
91  """Config for FitTanSipWcsTask."""
92  pass
93 
94 
95 class FitAffineWcsTask(pipeBase.Task):
96  """Fit a TAN-SIP WCS given a list of reference object/source matches.
97 
98  This WCS fitter should be used on top of a cameraGeom distortion model as
99  the model assumes that only a shift the WCS center position and a small
100  affine transform are required.
101  """
102  ConfigClass = FitAffineWcsConfig
103  _DefaultName = "fitAffineWcs"
104 
105  @timeMethod
106  def fitWcs(self,
107  matches,
108  initWcs,
109  bbox=None,
110  refCat=None,
111  sourceCat=None,
112  exposure=None):
113  """Fit a simple Affine transform with a shift to the matches and update
114  the WCS.
115 
116  This method assumes that the distortion model of the telescope is
117  applied correctly and is accurate with only a slight rotation,
118  rotation, and "squish" required to fit to the reference locations.
119 
120  Parameters
121  ----------
122  matches : `list` of `lsst.afw.table.ReferenceMatch`
123  The following fields are read:
124 
125  - match.first (reference object) coord
126  - match.second (source) centroid
127 
128  The following fields are written:
129 
130  - match.first (reference object) centroid,
131  - match.second (source) centroid
132  - match.distance (on sky separation, in radians)
133 
134  initWcs : `lsst.afw.geom.SkyWcs`
135  initial WCS
136  bbox : `lsst.geom.Box2I`
137  Ignored; present for consistency with FitSipDistortionTask.
138  refCat : `lsst.afw.table.SimpleCatalog`
139  reference object catalog, or None.
140  If provided then all centroids are updated with the new WCS,
141  otherwise only the centroids for ref objects in matches are
142  updated. Required fields are "centroid_x", "centroid_y",
143  "coord_ra", and "coord_dec".
144  sourceCat : `lsst.afw.table.SourceCatalog`
145  source catalog, or None.
146  If provided then coords are updated with the new WCS;
147  otherwise only the coords for sources in matches are updated.
148  Required fields are "slot_Centroid_x", "slot_Centroid_y", and
149  "coord_ra", and "coord_dec".
150  exposure : `lsst.afw.image.Exposure`
151  Ignored; present for consistency with FitSipDistortionTask.
152 
153  Returns
154  -------
155  result : `lsst.pipe.base.Struct`
156  with the following fields:
157 
158  - ``wcs`` : the fit WCS (`lsst.afw.geom.SkyWcs`)
159  - ``scatterOnSky`` : median on-sky separation between reference
160  objects and sources in "matches" (`lsst.afw.geom.Angle`)
161  """
162  # Create a data-structure that decomposes the input Wcs frames and
163  # appends the new transform.
164  wcsMaker = TransformedSkyWcsMaker(initWcs)
165 
166  refPoints = []
167  srcPixels = []
168  offsetLong = 0
169  offsetLat = 0
170  # Grab reference coordinates and source centroids. Compute the average
171  # direction and separation between the reference and the sources.
172  for match in matches:
173  refCoord = match.first.getCoord()
174  refPoints.append(refCoord)
175  srcCentroid = match.second.getCentroid()
176  srcPixels.append(srcCentroid)
177  srcCoord = initWcs.pixelToSky(srcCentroid)
178  deltaLong, deltaLat = srcCoord.getTangentPlaneOffset(refCoord)
179  offsetLong += deltaLong.asArcseconds()
180  offsetLat += deltaLat.asArcseconds()
181  offsetLong /= len(srcPixels)
182  offsetLat /= len(srcPixels)
183  offsetDist = np.sqrt(offsetLong ** 2 + offsetLat ** 2)
184  if offsetDist > 0.:
185  offsetDir = np.degrees(np.arccos(offsetLong / offsetDist))
186  else:
187  offsetDir = 0.
188  offsetDir *= np.sign(offsetLat)
189  self.log.debug("Initial shift guess: Direction: %.3f, Dist %.3f...",
190  offsetDir, offsetDist)
191 
192  # Best performing fitter in scipy tried so far (vs. default settings in
193  # minimize). Exits early because of the xTol value which cannot be
194  # disabled in scipy1.2.1. Matrix starting values are non-zero as this
195  # results in better fit off-diagonal terms.
196  fit = least_squares(
197  _chiFunc,
198  x0=[offsetDir, offsetDist, 1., 1e-8, 1e-8, 1.],
199  args=(refPoints, srcPixels, wcsMaker),
200  method='dogbox',
201  bounds=[[-360, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf],
202  [360, np.inf, np.inf, np.inf, np.inf, np.inf]],
203  ftol=2.3e-16,
204  gtol=2.31e-16,
205  xtol=2.3e-16)
206  self.log.debug("Best fit: Direction: %.3f, Dist: %.3f, "
207  "Affine matrix: [[%.6f, %.6f], [%.6f, %.6f]]...",
208  fit.x[0], fit.x[1],
209  fit.x[2], fit.x[3], fit.x[4], fit.x[5])
210 
211  wcs = wcsMaker.makeWcs(fit.x[:2], fit.x[2:].reshape((2, 2)))
212 
213  # Copied from other fit*WcsTasks.
214  if refCat is not None:
215  self.log.debug("Updating centroids in refCat")
216  lsst.afw.table.updateRefCentroids(wcs, refList=refCat)
217  else:
218  self.log.warning("Updating reference object centroids in match list; refCat is None")
220  wcs,
221  refList=[match.first for match in matches])
222 
223  if sourceCat is not None:
224  self.log.debug("Updating coords in sourceCat")
225  lsst.afw.table.updateSourceCoords(wcs, sourceList=sourceCat)
226  else:
227  self.log.warning("Updating source coords in match list; sourceCat is None")
229  wcs,
230  sourceList=[match.second for match in matches])
231  setMatchDistance(matches)
232 
233  stats = makeMatchStatisticsInRadians(wcs,
234  matches,
235  lsst.afw.math.MEDIAN)
236  scatterOnSky = stats.getValue() * radians
237 
238  self.log.debug("In fitter scatter %.4f", scatterOnSky.asArcseconds())
239 
240  return lsst.pipe.base.Struct(
241  wcs=wcs,
242  scatterOnSky=scatterOnSky,
243  )
244 
245 
247  """Convenience class for appending a shifting an input SkyWcs on sky and
248  appending an affine transform.
249 
250  The class assumes that all frames are sequential and are mapped one to the
251  next.
252 
253  Parameters
254  ----------
255  input_sky_wcs : `lsst.afw.geom.SkyWcs`
256  WCS to decompose and append affine matrix and shift in on sky
257  location to.
258  """
259 
260  def __init__(self, inputSkyWcs):
261  self.frameDictframeDict = inputSkyWcs.getFrameDict()
262 
263  # Grab the order of the frames by index.
264  # TODO: DM-20825
265  # Change the frame the transform is appended to to be explicitly
266  # the FIELD_ANGLE->IWC transform. Requires related tickets to be
267  # completed.
268  domains = self.frameDictframeDict.getAllDomains()
269  self.frameIdxsframeIdxs = np.sort([self.frameDictframeDict.getIndex(domain)
270  for domain in domains])
271  self.frameMinframeMin = np.min(self.frameIdxsframeIdxs)
272  self.frameMaxframeMax = np.max(self.frameIdxsframeIdxs)
273 
274  # Find frame just before the final mapping to sky and store those
275  # indices and mappings for later.
276  self.mapFrommapFrom = self.frameMaxframeMax - 2
277  if self.mapFrommapFrom < self.frameMinframeMin:
278  self.mapFrommapFrom = self.frameMinframeMin
279  self.mapTomapTo = self.frameMaxframeMax - 1
280  if self.mapTomapTo <= self.mapFrommapFrom:
281  self.mapTomapTo = self.frameMaxframeMax
282  self.lastMapBeforeSkylastMapBeforeSky = self.frameDictframeDict.getMapping(
283  self.mapFrommapFrom, self.mapTomapTo)
284 
285  # Get the original WCS sky location.
286 
287  self.originorigin = inputSkyWcs.getSkyOrigin()
288 
289  def makeWcs(self, crvalOffset, affMatrix):
290  """Apply a shift and affine transform to the WCS internal to this
291  class.
292 
293  A new SkyWcs with these transforms applied is returns.
294 
295  Parameters
296  ----------
297  crval_shift : `numpy.ndarray`, (2,)
298  Shift in radians to apply to the Wcs origin/crvals.
299  aff_matrix : 'numpy.ndarray', (3, 3)
300  Affine matrix to apply to the mapping/transform to add to the
301  WCS.
302 
303  Returns
304  -------
305  outputWcs : `lsst.afw.geom.SkyWcs`
306  Wcs with a final shift and affine transform applied.
307  """
308  # Create a WCS that only maps from IWC to Sky with the shifted
309  # Sky origin position. This is simply the final undistorted tangent
310  # plane to sky. The PIXELS to SKY map will be become our IWC to SKY
311  # map and gives us our final shift position.
312  iwcsToSkyWcs = makeSkyWcs(
313  Point2D(0., 0.),
314  self.originorigin.offset(crvalOffset[0] * degrees,
315  crvalOffset[1] * arcseconds),
316  np.array([[1., 0.], [0., 1.]]))
317  iwcToSkyMap = iwcsToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY")
318 
319  # Append a simple affine Matrix transform to the current to the
320  # second to last frame mapping. e.g. the one just before IWC to SKY.
321  newMapping = self.lastMapBeforeSkylastMapBeforeSky.then(astshim.MatrixMap(affMatrix))
322 
323  # Create a new frame dict starting from the input_sky_wcs's first
324  # frame. Append the correct mapping created above and our new on
325  # sky location.
326  outputFrameDict = astshim.FrameDict(
327  self.frameDictframeDict.getFrame(self.frameMinframeMin))
328  for frameIdx in self.frameIdxsframeIdxs:
329  if frameIdx == self.mapFrommapFrom:
330  outputFrameDict.addFrame(
331  self.mapFrommapFrom,
332  newMapping,
333  self.frameDictframeDict.getFrame(self.mapTomapTo))
334  elif frameIdx >= self.mapTomapTo:
335  continue
336  else:
337  outputFrameDict.addFrame(
338  frameIdx,
339  self.frameDictframeDict.getMapping(frameIdx, frameIdx + 1),
340  self.frameDictframeDict.getFrame(frameIdx + 1))
341  # Append the final sky frame to the frame dict.
342  outputFrameDict.addFrame(
343  self.frameMaxframeMax - 1,
344  iwcToSkyMap,
345  iwcsToSkyWcs.getFrameDict().getFrame("SKY"))
346 
347  return SkyWcs(outputFrameDict)
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None)
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
Point< double, 2 > Point2D
Definition: Point.h:324
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.