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
tractInfo.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 #
23 import lsst.afw.coord as afwCoord
24 import lsst.afw.geom as afwGeom
25 import lsst.afw.image as afwImage
26 from .patchInfo import PatchInfo
27 
28 __all__ = ["TractInfo"]
29 
30 class TractInfo(object):
31  """Information about a tract in a SkyMap sky pixelization
32 
33  The tract is subdivided into rectangular patches. Each patch has the following properties:
34  - An inner region defined by an inner bounding. The inner regions of the patches exactly tile the tract,
35  and all inner regions have the same dimensions. The tract is made larger as required to make this work.
36  - An outer region defined by an outer bounding box. The outer region extends beyond the inner region
37  by patchBorder pixels in all directions, except there is no border at the edges of the tract.
38  Thus patches overlap each other but never extend off the tract. If you do not want any overlap
39  between adjacent patches then set patchBorder to 0.
40  - An index that consists of a pair of integers:
41  0 <= x index < numPatches[0]
42  0 <= y index < numPatches[1]
43  Patch 0,0 is at the minimum corner of the tract bounding box.
44  """
45  def __init__(self, id, patchInnerDimensions, patchBorder, ctrCoord, vertexCoordList, tractOverlap, wcs):
46  """Construct a TractInfo
47 
48  @param[in] id: tract ID
49  @param[in] patchInnerDimensions: dimensions of inner region of patches (x,y pixels)
50  @param[in] patchBorder: overlap between adjacent patches (in pixels, one int)
51  @param[in] ctrCoord: sky coordinate of center of inner region of tract, as an afwCoord.Coord;
52  also used as the CRVAL for the WCS.
53  @param[in] vertexCoordList: list of sky coordinates (afwCoord.Coord)
54  of vertices that define the boundaries of the inner region
55  @param[in] tractOverlap: minimum overlap between adjacent sky tracts; an afwGeom.Angle;
56  this defines the minimum distance the tract extends beyond the inner region in all directions
57  @param[in,out] wcs: an afwImage.Wcs; the reference pixel will be shifted as required
58  so that the lower left-hand pixel (index 0,0) has pixel position 0.0, 0.0
59 
60  @warning
61  - It is not enforced that ctrCoord is the center of vertexCoordList, but SkyMap relies on it
62  - vertexCoordList will likely become a geom SphericalConvexPolygon someday.
63  """
64  self._id = id
65  try:
66  assert len(patchInnerDimensions) == 2
67  self._patchInnerDimensions = afwGeom.Extent2I(*(int(val) for val in patchInnerDimensions))
68  except:
69  raise TypeError("patchInnerDimensions=%s; must be two ints" % (patchInnerDimensions,))
70  self._patchBorder = int(patchBorder)
71  self._ctrCoord = ctrCoord
72  self._vertexCoordList = tuple(coord.clone() for coord in vertexCoordList)
73  self._tractOverlap = tractOverlap
74 
75  minBBox = self._minimumBoundingBox(wcs)
76  initialBBox, self._numPatches = self._setupPatches(minBBox, wcs)
77  self._bbox, self._wcs = self._finalOrientation(initialBBox, wcs)
78 
79 
80  def _minimumBoundingBox(self, wcs):
81  """Calculate the minimum bounding box for the tract, given the WCS
82 
83  The bounding box is created in the frame of the supplied WCS,
84  so that it's OK if the coordinates are negative.
85 
86  We compute the bounding box that holds all the vertices and the
87  desired overlap.
88  """
89  minBBoxD = afwGeom.Box2D()
90  halfOverlap = self._tractOverlap / 2.0
91  for vertexCoord in self._vertexCoordList:
92  vertexDeg = vertexCoord.getPosition(afwGeom.degrees)
93  if self._tractOverlap == 0:
94  minBBoxD.include(wcs.skyToPixel(vertexCoord))
95  else:
96  numAngles = 24
97  angleIncr = afwGeom.Angle(360.0, afwGeom.degrees) / float(numAngles)
98  for i in range(numAngles):
99  offAngle = angleIncr * i
100  offCoord = vertexCoord.clone()
101  offCoord.offset(offAngle, halfOverlap)
102  pixPos = wcs.skyToPixel(offCoord)
103  minBBoxD.include(pixPos)
104  return minBBoxD
105 
106  def _setupPatches(self, minBBox, wcs):
107  """Setup for patches of a particular size.
108 
109  We grow the bounding box to hold an exact multiple of
110  the desired size (patchInnerDimensions), while keeping
111  the center roughly the same. We return the final
112  bounding box, and the number of patches in each dimension
113  (as an Extent2I).
114 
115  @param minBBox Minimum bounding box for tract
116  @param wcs Wcs object
117  @return final bounding box, number of patches
118  """
119  bbox = afwGeom.Box2I(minBBox)
120  bboxMin = bbox.getMin()
121  bboxDim = bbox.getDimensions()
122  numPatches = afwGeom.Extent2I(0, 0)
123  for i, innerDim in enumerate(self._patchInnerDimensions):
124  num = (bboxDim[i] + innerDim - 1) // innerDim # round up
125  deltaDim = (innerDim * num) - bboxDim[i]
126  if deltaDim > 0:
127  bboxDim[i] = innerDim * num
128  bboxMin[i] -= deltaDim // 2
129  numPatches[i] = num
130  bbox = afwGeom.Box2I(bboxMin, bboxDim)
131  return bbox, numPatches
132 
133 
134  def _finalOrientation(self, bbox, wcs):
135  """Determine the final orientation
136 
137  We offset everything so the lower-left corner is at 0,0
138  and compute the final Wcs.
139 
140  @param bbox Current bounding box
141  @param wcs Current Wcs
142  @return revised bounding box, revised Wcs
143  """
144  finalBBox = afwGeom.Box2I(afwGeom.Point2I(0, 0), bbox.getDimensions())
145  # shift the WCS by the same amount as the bbox; extra code is required
146  # because simply subtracting makes an Extent2I
147  pixPosOffset = afwGeom.Extent2D(finalBBox.getMinX() - bbox.getMinX(),
148  finalBBox.getMinY() - bbox.getMinY())
149  wcs.shiftReferencePixel(pixPosOffset)
150  return finalBBox, wcs
151 
152 
153  def findPatch(self, coord):
154  """Find the patch containing the specified coord
155 
156  @param[in] coord: sky coordinate (afwCoord.Coord)
157  @return PatchInfo of patch whose inner bbox contains the specified coord
158 
159  @raise LookupError if coord is not in tract
160 
161  @note This routine will be more efficient if coord is ICRS.
162  """
163  pixelInd = afwGeom.Point2I(self.getWcs().skyToPixel(coord.toIcrs()))
164  if not self.getBBox().contains(pixelInd):
165  raise LookupError("coord %s is not in tract %s" % (coord, self.getId()))
166  patchInd = tuple(int(pixelInd[i]/self._patchInnerDimensions[i]) for i in range(2))
167  return self.getPatchInfo(patchInd)
168 
169  def findPatchList(self, coordList):
170  """Find patches containing the specified list of coords
171 
172  @param[in] coordList: list of sky coordinates (afwCoord.Coord)
173  @return list of PatchInfo for patches that contain, or may contain, the specified region.
174  The list will be empty if there is no overlap.
175 
176  @warning:
177  * This may give incorrect answers on regions that are larger than a tract
178  * This uses a naive algorithm that may find some patches that do not overlap the region
179  (especially if the region is not a rectangle aligned along patch x,y).
180  """
181  box2D = afwGeom.Box2D()
182  for coord in coordList:
183  try:
184  pixelPos = self.getWcs().skyToPixel(coord.toIcrs())
186  # the point is so far off the tract that its pixel position cannot be computed
187  continue
188  box2D.include(pixelPos)
189  bbox = afwGeom.Box2I(box2D)
190  bbox.grow(self.getPatchBorder())
191  bbox.clip(self.getBBox())
192  if bbox.isEmpty():
193  return ()
194 
195  llPatchInd = tuple(int(bbox.getMin()[i]/self._patchInnerDimensions[i]) for i in range(2))
196  urPatchInd = tuple(int(bbox.getMax()[i]/self._patchInnerDimensions[i]) for i in range(2))
197  return tuple(self.getPatchInfo((xInd, yInd))
198  for xInd in range(llPatchInd[0], urPatchInd[0]+1)
199  for yInd in range(llPatchInd[1], urPatchInd[1]+1))
200 
201  def getBBox(self):
202  """Get bounding box of tract (as an afwGeom.Box2I)
203  """
204  return afwGeom.Box2I(self._bbox)
205 
206  def getCtrCoord(self):
207  """Get sky coordinate of center of tract (as an afwCoord.Coord)
208  """
209  return self._ctrCoord
210 
211  def getId(self):
212  """Get ID of tract
213  """
214  return self._id
215 
216  def getNumPatches(self):
217  """Get the number of patches in x, y
218 
219  @return the number of patches in x, y
220  """
221  return self._numPatches
222 
223  def getPatchBorder(self):
224  """Get batch border
225 
226  @return patch border (pixels)
227  """
228  return self._patchBorder
229 
230  def getPatchInfo(self, index):
231  """Return information for the specified patch
232 
233  @param[in] index: index of patch, as a pair of ints
234  @return patch info, an instance of PatchInfo
235 
236  @raise IndexError if index is out of range
237  """
238  if (not 0 <= index[0] < self._numPatches[0]) \
239  or (not 0 <= index[1] < self._numPatches[1]):
240  raise IndexError("Patch index %s is not in range [0-%d, 0-%d]" % \
241  (index, self._numPatches[0]-1, self._numPatches[1]-1))
242  innerMin = afwGeom.Point2I(*[index[i] * self._patchInnerDimensions[i] for i in range(2)])
243  innerBBox = afwGeom.Box2I(innerMin, self._patchInnerDimensions)
244  if not self._bbox.contains(innerBBox):
245  raise RuntimeError(
246  "Bug: patch index %s valid but inner bbox=%s not contained in tract bbox=%s" % \
247  (index, innerBBox, self._bbox))
248  outerBBox = afwGeom.Box2I(innerBBox)
249  outerBBox.grow(self.getPatchBorder())
250  outerBBox.clip(self._bbox)
251  return PatchInfo(
252  index = index,
253  innerBBox = innerBBox,
254  outerBBox = outerBBox,
255  )
256 
258  """Get dimensions of inner region of the patches (all are the same)
259 
260  @return dimensions of inner region of the patches (as an afwGeom Extent2I)
261  """
262  return self._patchInnerDimensions
263 
264  def getTractOverlap(self):
265  """Get minimum overlap of adjacent sky tracts
266 
267  @return minimum overlap between adjacent sky tracts, as an afwGeom Angle
268  """
269  return self._tractOverlap
270 
271  def getVertexList(self):
272  """Get list of sky coordinates of vertices that define the boundary of the inner region
273 
274  @warning: this is not a deep copy
275  @warning vertexCoordList will likely become a geom SphericalConvexPolygon someday.
276  """
277  return self._vertexCoordList
278 
279  def getWcs(self):
280  """Get WCS of tract
281 
282  @warning: this is not a deep copy
283  """
284  return self._wcs
285 
286  def __str__(self):
287  return "TractInfo(id=%s)" % (self._id,)
288 
289  def __repr__(self):
290  return "TractInfo(id=%s, ctrCoord=%s)" % (self._id, self._ctrCoord.getVector())
291 
292  def __iter__(self):
293  xNum, yNum = self.getNumPatches()
294  for y in range(yNum):
295  for x in range(xNum):
296  yield self.getPatchInfo((x,y))
297 
298  def __len__(self):
299  xNum, yNum = self.getNumPatches()
300  return xNum*yNum
301 
302  def __getitem__(self, index):
303  return self.getPatchInfo(index)
304 
305 
307  """Information for a tract specified explicitly
308 
309  A tract is placed at the explicitly defined coordinates, with the nominated
310  radius. The tracts are square (i.e., the radius is really a half-size).
311  """
312  def __init__(self, ident, patchInnerDimensions, patchBorder, ctrCoord, radius, tractOverlap, wcs):
313  vertexList = []
314  self._radius = radius
315  super(ExplicitTractInfo, self).__init__(ident, patchInnerDimensions, patchBorder, ctrCoord,
316  vertexList, tractOverlap, wcs)
317 
318  def _minimumBoundingBox(self, wcs):
319  """The minimum bounding box is calculated using the nominated radius"""
320  bbox = afwGeom.Box2D()
321  for i in range(4):
322  coord = self._ctrCoord.clone()
323  coord.offset(i * 90 * afwGeom.degrees, self._radius + self._tractOverlap)
324  pixPos = wcs.skyToPixel(coord)
325  bbox.include(pixPos)
326  return bbox
327 
An integer coordinate rectangle.
Definition: Box.h:53
A floating-point coordinate rectangle geometry.
Definition: Box.h:271