Compare two WCS over a rectangular grid of pixel positions.
64 maxDiffPix=0.01, nx=5, ny=5, doShortCircuit=
True):
65 """!Compare two WCS over a rectangular grid of pixel positions
67 @param[in] wcs0 WCS 0 (an lsst.afw.image.Wcs)
68 @param[in] wcs1 WCS 1 (an lsst.afw.image.Wcs)
69 @param[in] bbox boundaries of pixel grid over which to compare the WCSs (an lsst.afw.geom.Box2I or Box2D)
70 @param[in] maxDiffSky maximum separation between sky positions computed using Wcs.pixelToSky
71 (an lsst.afw.geom.Angle)
72 @param[in] maxDiffPix maximum separation between pixel positions computed using Wcs.skyToPixel
73 @param[in] nx number of points in x for the grid of pixel positions
74 @param[in] ny number of points in y for the grid of pixel positions
75 @param[in] doShortCircuit if True then stop at the first error, else test all values in the grid
76 and return information about the worst violations found
78 @return return an empty string if the WCS are sufficiently close; else return a string describing
79 the largest error measured in pixel coordinates (if sky to pixel error was excessive) and sky coordinates
80 (if pixel to sky error was excessive). If doShortCircuit is true then the reported error is likely to be
81 much less than the maximum error across the whole pixel grid.
84 raise RuntimeError(
"nx = %s and ny = %s must both be positive" % (nx, ny))
85 if maxDiffSky <= 0*afwGeom.arcseconds:
86 raise RuntimeError(
"maxDiffSky = %s must be positive" % (maxDiffSky,))
88 raise RuntimeError(
"maxDiffPix = %s must be positive" % (maxDiffPix,))
91 xList = numpy.linspace(bboxd.getMinX(), bboxd.getMaxX(), nx)
92 yList = numpy.linspace(bboxd.getMinY(), bboxd.getMaxY(), ny)
94 measDiffSky = (maxDiffSky,
"?")
95 measDiffPix = (maxDiffPix,
"?")
96 for x, y
in itertools.product(xList, yList):
98 sky0 = wcs0.pixelToSky(fromPixPos)
99 sky1 = wcs1.pixelToSky(fromPixPos)
100 diffSky = sky0.angularSeparation(sky1)
101 if diffSky > measDiffSky[0]:
102 measDiffSky = (diffSky, fromPixPos)
106 toPixPos0 = wcs0.skyToPixel(sky0)
107 toPixPos1 = wcs1.skyToPixel(sky0)
108 diffPix = math.hypot(*(toPixPos0 - toPixPos1))
109 if diffPix > measDiffPix[0]:
110 measDiffPix = (diffPix, sky0)
115 if measDiffSky[0] > maxDiffSky:
116 msgList.append(
"%s arcsec max measured sky error > %s arcsec max allowed sky error at pix pos=%s" %
117 (measDiffSky[0].asArcseconds(), maxDiffSky.asArcseconds(), measDiffSky[1]))
118 if measDiffPix[0] > maxDiffPix:
119 msgList.append(
"%s max measured pix error > %s max allowed pix error at sky pos=%s" %
120 (measDiffPix[0], maxDiffPix, measDiffPix[1]))
122 return "; ".join(msgList)