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
basicUtils.py
Go to the documentation of this file.
1 from __future__ import absolute_import, division
2 #
3 # LSST Data Management System
4 # Copyright 2008, 2009, 2010 LSST Corporation.
5 #
6 # This product includes software developed by the
7 # LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20 # the GNU General Public License along with this program. If not,
21 # see <http://www.lsstcorp.org/LegalNotices/>.
22 #
23 
24 """Application Framework image-related classes including Image, Mask and MaskedImage
25 """
26 import itertools
27 import math
28 
29 import numpy
30 
31 import lsst.utils.tests
32 import lsst.afw.geom as afwGeom
33 from . import imageLib
34 
35 __all__ = ["makeImageFromArray", "makeMaskFromArray", "makeMaskedImageFromArrays",
36  "wcsNearlyEqualOverBBox", "assertWcsNearlyEqualOverBBox"]
37 
38 suffixes = {str(numpy.uint16): "U", str(numpy.int32): "I", str(numpy.float32): "F", str(numpy.float64): "D"}
39 
40 def makeImageFromArray(array):
41  """Construct an Image from a NumPy array, inferring the Image type from the NumPy type.
42  Return None if input is None.
43  """
44  if array is None: return None
45  cls = getattr(imageLib, "Image%s" % (suffixes[str(array.dtype.type)],))
46  return cls(array)
47 
48 def makeMaskFromArray(array):
49  """Construct an Mask from a NumPy array, inferring the Mask type from the NumPy type.
50  Return None if input is None.
51  """
52  if array is None: return None
53  cls = getattr(imageLib, "Mask%s" % (suffixes[str(array.dtype.type)],))
54  return cls(array)
55 
56 def makeMaskedImageFromArrays(image, mask=None, variance=None):
57  """Construct a MaskedImage from three NumPy arrays, inferring the MaskedImage types from the NumPy types.
58  """
59  cls = getattr(imageLib, "MaskedImage%s" % (suffixes[str(image.dtype.type)],))
60  return cls(makeImageFromArray(image), makeMaskFromArray(mask), makeImageFromArray(variance))
61 
62 def _compareWcsOverBBox(wcs0, wcs1, bbox, maxDiffSky=0.01*afwGeom.arcseconds,
63  maxDiffPix=0.01, nx=5, ny=5, doShortCircuit=True):
64  """!Compare two WCS over a rectangular grid of pixel positions
65 
66  @param[in] wcs0 WCS 0 (an lsst.afw.image.Wcs)
67  @param[in] wcs1 WCS 1 (an lsst.afw.image.Wcs)
68  @param[in] bbox boundaries of pixel grid over which to compare the WCSs (an lsst.afw.geom.Box2I or Box2D)
69  @param[in] maxDiffSky maximum separation between sky positions computed using Wcs.pixelToSky
70  (an lsst.afw.geom.Angle)
71  @param[in] maxDiffPix maximum separation between pixel positions computed using Wcs.skyToPixel
72  @param[in] nx number of points in x for the grid of pixel positions
73  @param[in] ny number of points in y for the grid of pixel positions
74  @param[in] doShortCircuit if True then stop at the first error, else test all values in the grid
75  and return information about the worst violations found
76 
77  @return return an empty string if the WCS are sufficiently close; else return a string describing
78  the largest error measured in pixel coordinates (if sky to pixel error was excessive) and sky coordinates
79  (if pixel to sky error was excessive). If doShortCircuit is true then the reported error is likely to be
80  much less than the maximum error across the whole pixel grid.
81  """
82  if nx < 1 or ny < 1:
83  raise RuntimeError("nx = %s and ny = %s must both be positive" % (nx, ny))
84  if maxDiffSky <= 0*afwGeom.arcseconds:
85  raise RuntimeError("maxDiffSky = %s must be positive" % (maxDiffSky,))
86  if maxDiffPix <= 0:
87  raise RuntimeError("maxDiffPix = %s must be positive" % (maxDiffPix,))
88 
89  bboxd = afwGeom.Box2D(bbox)
90  xList = numpy.linspace(bboxd.getMinX(), bboxd.getMaxX(), nx)
91  yList = numpy.linspace(bboxd.getMinY(), bboxd.getMaxY(), ny)
92  # we don't care about measured error unless it is too large, so initialize to max allowed
93  measDiffSky = (maxDiffSky, "?") # (sky diff, pix pos)
94  measDiffPix = (maxDiffPix, "?") # (pix diff, sky pos)
95  for x, y in itertools.product(xList, yList):
96  fromPixPos = afwGeom.Point2D(x, y)
97  sky0 = wcs0.pixelToSky(fromPixPos)
98  sky1 = wcs1.pixelToSky(fromPixPos)
99  diffSky = sky0.angularSeparation(sky1)
100  if diffSky > measDiffSky[0]:
101  measDiffSky = (diffSky, fromPixPos)
102  if doShortCircuit:
103  break
104 
105  toPixPos0 = wcs0.skyToPixel(sky0)
106  toPixPos1 = wcs1.skyToPixel(sky0)
107  diffPix = math.hypot(*(toPixPos0 - toPixPos1))
108  if diffPix > measDiffPix[0]:
109  measDiffPix = (diffPix, sky0)
110  if doShortCircuit:
111  break
112 
113  msgList = []
114  if measDiffSky[0] > maxDiffSky:
115  msgList.append("%s arcsec max measured sky error > %s arcsec max allowed sky error at pix pos=%s" %
116  (measDiffSky[0].asArcseconds(), maxDiffSky.asArcseconds(), measDiffSky[1]))
117  if measDiffPix[0] > maxDiffPix:
118  msgList.append("%s max measured pix error > %s max allowed pix error at sky pos=%s" %
119  (measDiffPix[0], maxDiffPix, measDiffPix[1]))
120 
121  return "; ".join(msgList)
122 
123 def wcsNearlyEqualOverBBox(wcs0, wcs1, bbox, maxDiffSky=0.01*afwGeom.arcseconds,
124  maxDiffPix=0.01, nx=5, ny=5):
125  """!Return True if two WCS are nearly equal over a grid of pixel positions, else False
126 
127  @param[in] wcs0 WCS 0 (an lsst.afw.image.Wcs)
128  @param[in] wcs1 WCS 1 (an lsst.afw.image.Wcs)
129  @param[in] bbox boundaries of pixel grid over which to compare the WCSs (an lsst.afw.geom.Box2I or Box2D)
130  @param[in] maxDiffSky maximum separation between sky positions computed using Wcs.pixelToSky
131  (an lsst.afw.geom.Angle)
132  @param[in] maxDiffPix maximum separation between pixel positions computed using Wcs.skyToPixel
133  @param[in] nx number of points in x for the grid of pixel positions
134  @param[in] ny number of points in y for the grid of pixel positions
135  @param[in] doShortCircuit if True then stop at the first error, else test all values in the grid
136  and return information about the worst violations found
137  """
138  return not bool(_compareWcsOverBBox(
139  wcs0 = wcs0,
140  wcs1 = wcs1,
141  bbox = bbox,
142  maxDiffSky = maxDiffSky,
143  maxDiffPix = maxDiffPix,
144  nx = nx,
145  ny = ny,
146  doShortCircuit = True,
147  ))
148 
149 @lsst.utils.tests.inTestCase
150 def assertWcsNearlyEqualOverBBox(testCase, wcs0, wcs1, bbox, maxDiffSky=0.01*afwGeom.arcseconds, maxDiffPix=0.01,
151  nx=5, ny=5, msg="WCSs differ"):
152  """!Compare pixelToSky and skyToPixel for two WCS over a rectangular grid of pixel positions
153 
154  If the WCS are too divergent, call testCase.fail; the message describes the largest error measured
155  in pixel coordinates (if sky to pixel error was excessive) and sky coordinates (if pixel to sky error
156  was excessive) across the entire pixel grid.
157 
158  @param[in] testCase unittest.TestCase instance the test is part of;
159  an object supporting one method: fail(self, msgStr)
160  @param[in] wcs0 WCS 0 (an lsst.afw.image.Wcs)
161  @param[in] wcs1 WCS 1 (an lsst.afw.image.Wcs)
162  @param[in] bbox boundaries of pixel grid over which to compare the WCSs (an lsst.afw.geom.Box2I or Box2D)
163  @param[in] maxDiffSky maximum separation between sky positions computed using Wcs.pixelToSky
164  (an lsst.afw.geom.Angle)
165  @param[in] maxDiffPix maximum separation between pixel positions computed using Wcs.skyToPixel
166  @param[in] nx number of points in x for the grid of pixel positions
167  @param[in] ny number of points in y for the grid of pixel positions
168  @param[in] msg exception message prefix; details of the error are appended after ": "
169  """
170  errMsg = _compareWcsOverBBox(
171  wcs0 = wcs0,
172  wcs1 = wcs1,
173  bbox = bbox,
174  maxDiffSky = maxDiffSky,
175  maxDiffPix = maxDiffPix,
176  nx = nx,
177  ny = ny,
178  doShortCircuit = False,
179  )
180  if errMsg:
181  testCase.fail("%s: %s" % (msg, errMsg))
def wcsNearlyEqualOverBBox
Return True if two WCS are nearly equal over a grid of pixel positions, else False.
Definition: basicUtils.py:124
def _compareWcsOverBBox
Compare two WCS over a rectangular grid of pixel positions.
Definition: basicUtils.py:63
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