39 maxDiffPix=0.01, nx=5, ny=5, doShortCircuit=True):
40 """Compare two :py:class:`WCS <lsst.afw.geom.SkyWcs>` over a rectangular grid of pixel positions
41
42 Parameters
43 ----------
44 wcs0 : `lsst.afw.geom.SkyWcs`
45 WCS 0
46 wcs1 : `lsst.afw.geom.SkyWcs`
47 WCS 1
48 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D`
49 boundaries of pixel grid over which to compare the WCSs
50 maxDiffSky : `lsst.geom.Angle`
51 maximum separation between sky positions computed using Wcs.pixelToSky
52 maxDiffPix : `float`
53 maximum separation between pixel positions computed using Wcs.skyToPixel
54 nx : `int`
55 number of points in x for the grid of pixel positions
56 ny : `int`
57 number of points in y for the grid of pixel positions
58 doShortCircuit : `bool`
59 if True then stop at the first error, else test all values in the grid
60 and return information about the worst violations found
61
62 Returns
63 -------
64 msg : `str`
65 an empty string if the WCS are sufficiently close; else return a string describing
66 the largest error measured in pixel coordinates (if sky to pixel error was excessive)
67 and sky coordinates (if pixel to sky error was excessive). If doShortCircuit is true
68 then the reported error is likely to be much less than the maximum error across the
69 whole pixel grid.
70 """
71 if nx < 1 or ny < 1:
72 raise RuntimeError(f"nx = {nx} and ny = {ny} must both be positive")
73 if maxDiffSky < 0*lsst.geom.arcseconds:
74 raise RuntimeError(f"maxDiffSky = {maxDiffSky} must not be negative")
75 if maxDiffPix < 0:
76 raise RuntimeError(f"maxDiffPix = {maxDiffPix} must not be negative")
77
79 xList = np.linspace(bboxd.getMinX(), bboxd.getMaxX(), nx)
80 yList = np.linspace(bboxd.getMinY(), bboxd.getMaxY(), ny)
81
82
83 measDiffSky = (maxDiffSky, "?")
84 measDiffPix = (maxDiffPix, "?")
85 for x, y in itertools.product(xList, yList):
87 sky0 = wcs0.pixelToSky(fromPixPos)
88 sky1 = wcs1.pixelToSky(fromPixPos)
89 diffSky = sky0.separation(sky1)
90 if diffSky > measDiffSky[0]:
91 measDiffSky = (diffSky, fromPixPos)
92 if doShortCircuit:
93 break
94
95 toPixPos0 = wcs0.skyToPixel(sky0)
96 toPixPos1 = wcs1.skyToPixel(sky0)
97 diffPix = math.hypot(*(toPixPos0 - toPixPos1))
98 if diffPix > measDiffPix[0]:
99 measDiffPix = (diffPix, sky0)
100 if doShortCircuit:
101 break
102
103 msgList = []
104 if measDiffSky[0] > maxDiffSky:
105 msgList.append(f"{measDiffSky[0].asArcseconds()} arcsec max measured sky error "
106 f"> {maxDiffSky.asArcseconds()} arcsec max allowed sky error "
107 f"at pix pos=measDiffSky[1]")
108 if measDiffPix[0] > maxDiffPix:
109 msgList.append(f"{measDiffPix[0]} max measured pix error "
110 f"> {maxDiffPix} max allowed pix error "
111 f"at sky pos={measDiffPix[1]}")
112
113 return "; ".join(msgList)
114
115
A floating-point coordinate rectangle geometry.