LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
approximateWcs.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 
24 import numpy
25 import lsst.afw.image as afwImage
26 import lsst.afw.table as afwTable
27 import lsst.afw.geom as afwGeom
28 from lsst.meas.base import SingleFrameMeasurementTask
29 from lsst.meas.astrom.sip import makeCreateWcsWithSip
30 from lsst.afw.image.basicUtils import assertWcsNearlyEqualOverBBox
31 
32 __all__ = ["approximateWcs"]
33 
34 class _MockTestCase(object):
35  """A fake unit test case class that will enable us to call
36  assertWcsNearlyEqualOverBBox from the method approximateWcs"""
37 
38  def fail(self, msgStr):
39  raise UserWarning("WCS fitting failed " + msgStr)
40 
41 def approximateWcs(wcs, bbox, order=3, nx=20, ny=20, iterations=3, \
42  skyTolerance=0.001*afwGeom.arcseconds, pixelTolerance=0.02, useTanWcs=False):
43  """Approximate an existing WCS as a TAN-SIP WCS
44 
45  The fit is performed by evaluating the WCS at a uniform grid of points within a bounding box.
46 
47  @param[in] wcs wcs to approximate
48  @param[in] bbox the region over which the WCS will be fit
49  @param[in] order order of SIP fit
50  @param[in] nx number of grid points along x
51  @param[in] ny number of grid points along y
52  @param[in] iterations number of times to iterate over fitting
53  @param[in] skyTolerance maximum allowed difference in world coordinates between
54  input wcs and approximate wcs (default is 0.001 arcsec)
55  @param[in] pixelTolerance maximum allowed difference in pixel coordinates between
56  input wcs and approximate wcs (default is 0.02 pixels)
57  @param[in] useTanWcs send a TAN version of wcs to the fitter? It is documented to require that,
58  but I don't think the fitter actually cares
59  @return the fit TAN-SIP WCS
60  """
61  if useTanWcs:
62  crCoord = wcs.getSkyOrigin()
63  crPix = wcs.getPixelOrigin()
64  cdMat = wcs.getCDMatrix()
65  tanWcs = afwImage.makeWcs(crCoord, crPix, cdMat[0,0], cdMat[0,1], cdMat[1,0], cdMat[1,1])
66  else:
67  tanWcs = wcs
68 
69  # create a matchList consisting of a grid of points covering the bbox
70  refSchema = afwTable.SimpleTable.makeMinimalSchema()
71  refCoordKey = afwTable.CoordKey(refSchema["coord"])
72  refCat = afwTable.SimpleCatalog(refSchema)
73 
74  sourceSchema = afwTable.SourceTable.makeMinimalSchema()
75  SingleFrameMeasurementTask(schema=sourceSchema) # expand the schema
76  sourceCentroidKey = afwTable.Point2DKey(sourceSchema["slot_Centroid"])
77 
78  sourceCat = afwTable.SourceCatalog(sourceSchema)
79 
80  matchList = afwTable.ReferenceMatchVector()
81 
82  bboxd = afwGeom.Box2D(bbox)
83  for x in numpy.linspace(bboxd.getMinX(), bboxd.getMaxX(), nx):
84  for y in numpy.linspace(bboxd.getMinY(), bboxd.getMaxY(), ny):
85  pixelPos = afwGeom.Point2D(x, y)
86  skyCoord = wcs.pixelToSky(pixelPos)
87 
88  refObj = refCat.addNew()
89  refObj.set(refCoordKey, skyCoord)
90 
91  source = sourceCat.addNew()
92  source.set(sourceCentroidKey, pixelPos)
93 
94  matchList.append(afwTable.ReferenceMatch(refObj, source, 0.0))
95 
96  # The TAN-SIP fitter is fitting x and y separately, so we have to iterate to make it converge
97  for indx in range(iterations) :
98  sipObject = makeCreateWcsWithSip(matchList, tanWcs, order, bbox)
99  tanWcs = sipObject.getNewWcs()
100  fitWcs = sipObject.getNewWcs()
101 
102  mockTest = _MockTestCase()
103  assertWcsNearlyEqualOverBBox(mockTest, wcs, fitWcs, bbox, maxDiffSky=skyTolerance,
104  maxDiffPix=pixelTolerance)
105 
106  return fitWcs
std::vector< ReferenceMatch > ReferenceMatchVector
Definition: fwd.h:97
CreateWcsWithSip< MatchT > makeCreateWcsWithSip(std::vector< MatchT > const &matches, afw::image::Wcs const &linearWcs, int const order, afw::geom::Box2I const &bbox=afw::geom::Box2I(), int const ngrid=0)
Factory function for CreateWcsWithSip.
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
Lightweight representation of a geometric match between two records.
Definition: fwd.h:90
Wcs::Ptr makeWcs(coord::Coord const &crval, geom::Point2D const &crpix, double CD11, double CD12, double CD21, double CD22)
Create a Wcs object from crval, crpix, CD, using CD elements (useful from python) ...
Definition: makeWcs.cc:141
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
def assertWcsNearlyEqualOverBBox
Compare pixelToSky and skyToPixel for two WCS over a rectangular grid of pixel positions.
Definition: basicUtils.py:151
A FunctorKey used to get or set celestial coordiantes from a pair of Angle keys.
Definition: aggregates.h:119