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