LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
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 #
22 
23 __all__ = ["TractInfo"]
24 
25 import numbers
26 
28 import lsst.geom as geom
29 from lsst.sphgeom import ConvexPolygon
30 
31 from .patchInfo import PatchInfo, makeSkyPolygonFromBBox
32 
33 
34 class TractInfo:
35  """Information about a tract in a SkyMap sky pixelization
36 
37  Parameters
38  ----------
39  id : `int`
40  tract ID
41  patchInnerDimensions : `tuple` of `int`
42  Dimensions of inner region of patches (x,y pixels).
43  patchBorder : `int`
44  Overlap between adjacent patches (in pixels)
45  ctrCoord : `lsst.geom.SpherePoint`
46  ICRS sky coordinate of center of inner region of tract; also used as
47  the CRVAL for the WCS.
48  vertexCoordList : `list` of `lsst.geom.SpherePoint`
49  Vertices that define the boundaries of the inner region.
50  tractOverlap : `lsst.geom.Angle`
51  Minimum overlap between adjacent sky tracts; this defines the minimum
52  distance the tract extends beyond the inner region in all directions.
53  wcs : `lsst.afw.image.SkyWcs`
54  WCS for tract. The reference pixel will be shifted as required so that
55  the lower left-hand pixel (index 0,0) has pixel position 0.0, 0.0.
56 
57  Notes
58  -----
59  The tract is subdivided into rectangular patches. Each patch has the
60  following properties:
61 
62  - An inner region defined by an inner bounding box. The inner regions of
63  the patches exactly tile the tract, and all inner regions have the same
64  dimensions. The tract is made larger as required to make this work.
65 
66  - An outer region defined by an outer bounding box. The outer region
67  extends beyond the inner region by patchBorder pixels in all directions,
68  except there is no border at the edges of the tract.
69  Thus patches overlap each other but never extend off the tract.
70  If you do not want any overlap between adjacent patches then set
71  patchBorder to 0.
72 
73  - An index that consists of a pair of integers:
74 
75  * 0 <= x index < numPatches[0]
76 
77  * 0 <= y index < numPatches[1]
78 
79  Patch 0,0 is at the minimum corner of the tract bounding box.
80 
81  - It is not enforced that ctrCoord is the center of vertexCoordList, but
82  SkyMap relies on it.
83  """
84 
85  def __init__(self, id, patchInnerDimensions, patchBorder, ctrCoord, vertexCoordList, tractOverlap, wcs):
86  self._id_id = id
87  try:
88  assert len(patchInnerDimensions) == 2
89  self._patchInnerDimensions_patchInnerDimensions = geom.Extent2I(*(int(val) for val in patchInnerDimensions))
90  except Exception:
91  raise TypeError("patchInnerDimensions=%s; must be two ints" % (patchInnerDimensions,))
92  self._patchBorder_patchBorder = int(patchBorder)
93  self._ctrCoord_ctrCoord = ctrCoord
94  self._vertexCoordList_vertexCoordList = tuple(vertexCoordList)
95  self._tractOverlap_tractOverlap = tractOverlap
96 
97  minBBox = self._minimumBoundingBox_minimumBoundingBox(wcs)
98  initialBBox, self._numPatches_numPatches = self._setupPatches_setupPatches(minBBox, wcs)
99  self._bbox, self._wcs_wcs = self._finalOrientation_finalOrientation(initialBBox, wcs)
100 
101  def _minimumBoundingBox(self, wcs):
102  """Calculate the minimum bounding box for the tract, given the WCS.
103 
104  The bounding box is created in the frame of the supplied WCS,
105  so that it's OK if the coordinates are negative.
106 
107  We compute the bounding box that holds all the vertices and the
108  desired overlap.
109  """
110  minBBoxD = geom.Box2D()
111  halfOverlap = self._tractOverlap_tractOverlap / 2.0
112  for vertexCoord in self._vertexCoordList_vertexCoordList:
113  if self._tractOverlap_tractOverlap == 0:
114  minBBoxD.include(wcs.skyToPixel(vertexCoord))
115  else:
116  numAngles = 24
117  angleIncr = geom.Angle(360.0, geom.degrees) / float(numAngles)
118  for i in range(numAngles):
119  offAngle = angleIncr * i
120  offCoord = vertexCoord.offset(offAngle, halfOverlap)
121  pixPos = wcs.skyToPixel(offCoord)
122  minBBoxD.include(pixPos)
123  return minBBoxD
124 
125  def _setupPatches(self, minBBox, wcs):
126  """Setup for patches of a particular size.
127 
128  We grow the bounding box to hold an exact multiple of
129  the desired size (patchInnerDimensions), while keeping
130  the center roughly the same. We return the final
131  bounding box, and the number of patches in each dimension
132  (as an Extent2I).
133 
134  Parameters
135  ----------
136  minBBox : `lsst.geom.Box2I`
137  Minimum bounding box for tract
138  wcs : `lsst.afw.geom.SkyWcs`
139  Wcs object
140 
141  Returns
142  -------
143  bbox : `lsst.geom.Box2I
144  final bounding box, number of patches
145  numPatches : `int`
146  """
147  bbox = geom.Box2I(minBBox)
148  bboxMin = bbox.getMin()
149  bboxDim = bbox.getDimensions()
150  numPatches = geom.Extent2I(0, 0)
151  for i, innerDim in enumerate(self._patchInnerDimensions_patchInnerDimensions):
152  num = (bboxDim[i] + innerDim - 1) // innerDim # round up
153  deltaDim = (innerDim * num) - bboxDim[i]
154  if deltaDim > 0:
155  bboxDim[i] = innerDim * num
156  bboxMin[i] -= deltaDim // 2
157  numPatches[i] = num
158  bbox = geom.Box2I(bboxMin, bboxDim)
159  return bbox, numPatches
160 
161  def _finalOrientation(self, bbox, wcs):
162  """Determine the final orientation
163 
164  We offset everything so the lower-left corner is at 0,0
165  and compute the final Wcs.
166 
167  Parameters
168  ----------
169  bbox : `lsst.geom.Box2I`
170  Current bounding box.
171  wcs : `lsst.afw.geom.SkyWcs
172  Current Wcs.
173 
174  Returns
175  -------
176  finalBBox : `lsst.geom.Box2I`
177  Revised bounding box.
178  wcs : `lsst.afw.geom.SkyWcs`
179  Revised Wcs.
180  """
181  finalBBox = geom.Box2I(geom.Point2I(0, 0), bbox.getDimensions())
182  # shift the WCS by the same amount as the bbox; extra code is required
183  # because simply subtracting makes an Extent2I
184  pixPosOffset = geom.Extent2D(finalBBox.getMinX() - bbox.getMinX(),
185  finalBBox.getMinY() - bbox.getMinY())
186  wcs = wcs.copyAtShiftedPixelOrigin(pixPosOffset)
187  return finalBBox, wcs
188 
189  def getSequentialPatchIndex(self, patchInfo):
190  """Return a single integer that uniquely identifies the given patch
191  within this tract.
192  """
193  x, y = patchInfo.getIndex()
194  nx, ny = self.getNumPatchesgetNumPatches()
195  return nx*y + x
196 
197  def getPatchIndexPair(self, sequentialIndex):
198  nx, ny = self.getNumPatchesgetNumPatches()
199  x = sequentialIndex % nx
200  y = (sequentialIndex - x) / nx
201  return (x, y)
202 
203  def findPatch(self, coord):
204  """Find the patch containing the specified coord.
205 
206  Parameters
207  ----------
208  coord : `lsst.geom.SpherePoint`
209  ICRS sky coordinate to search for.
210 
211  Returns
212  -------
213  result : `lsst.skymap.PatchInfo`
214  PatchInfo of patch whose inner bbox contains the specified coord
215 
216  Raises
217  ------
218  LookupError
219  If coord is not in tract or we cannot determine the
220  pixel coordinate (which likely means the coord is off the tract).
221  """
222  try:
223  pixel = self.getWcsgetWcs().skyToPixel(coord)
225  # Point must be way off the tract
226  raise LookupError("Unable to determine pixel position for coordinate %s" % (coord,))
227  pixelInd = geom.Point2I(pixel)
228  if not self.getBBoxgetBBox().contains(pixelInd):
229  raise LookupError("coord %s is not in tract %s" % (coord, self.getIdgetId()))
230  patchInd = tuple(int(pixelInd[i]/self._patchInnerDimensions_patchInnerDimensions[i]) for i in range(2))
231  return self.getPatchInfogetPatchInfo(patchInd)
232 
233  def findPatchList(self, coordList):
234  """Find patches containing the specified list of coords.
235 
236  Parameters
237  ----------
238  coordList : `list` of `lsst.geom.SpherePoint`
239  ICRS sky coordinates to search for.
240 
241  Returns
242  -------
243  result : `list` of `lsst.skymap.PatchInfo`
244  List of PatchInfo for patches that contain, or may contain, the
245  specified region. The list will be empty if there is no overlap.
246 
247  Notes
248  -----
249  **Warning:**
250 
251  - This may give incorrect answers on regions that are larger than a
252  tract.
253 
254  - This uses a naive algorithm that may find some patches that do not
255  overlap the region (especially if the region is not a rectangle
256  aligned along patch x,y).
257  """
258  box2D = geom.Box2D()
259  for coord in coordList:
260  try:
261  pixelPos = self.getWcsgetWcs().skyToPixel(coord)
263  # the point is so far off the tract that its pixel position cannot be computed
264  continue
265  box2D.include(pixelPos)
266  bbox = geom.Box2I(box2D)
267  bbox.grow(self.getPatchBordergetPatchBorder())
268  bbox.clip(self.getBBoxgetBBox())
269  if bbox.isEmpty():
270  return ()
271 
272  llPatchInd = tuple(int(bbox.getMin()[i]/self._patchInnerDimensions_patchInnerDimensions[i]) for i in range(2))
273  urPatchInd = tuple(int(bbox.getMax()[i]/self._patchInnerDimensions_patchInnerDimensions[i]) for i in range(2))
274  return tuple(self.getPatchInfogetPatchInfo((xInd, yInd))
275  for xInd in range(llPatchInd[0], urPatchInd[0]+1)
276  for yInd in range(llPatchInd[1], urPatchInd[1]+1))
277 
278  def getBBox(self):
279  """Get bounding box of tract (as an geom.Box2I)
280  """
281  return geom.Box2I(self._bbox)
282 
283  def getCtrCoord(self):
284  """Get ICRS sky coordinate of center of tract
285  (as an lsst.geom.SpherePoint)
286  """
287  return self._ctrCoord_ctrCoord
288 
289  def getId(self):
290  """Get ID of tract
291  """
292  return self._id_id
293 
294  def getNumPatches(self):
295  """Get the number of patches in x, y.
296 
297  Returns
298  -------
299  result : `tuple` of `int`
300  The number of patches in x, y
301  """
302  return self._numPatches_numPatches
303 
304  def getPatchBorder(self):
305  return self._patchBorder_patchBorder
306 
307  def getPatchInfo(self, index):
308  """Return information for the specified patch.
309 
310  Parameters
311  ----------
312  index : `tuple` of `int`
313  Index of patch, as a pair of ints;
314  or a sequential index as returned by getSequentialPatchIndex;
315  negative values are not supported.
316 
317  Returns
318  -------
319  result : `lsst.skymap.PatchInfo`
320  The patch info for that index.
321 
322  Raises
323  ------
324  IndexError
325  If index is out of range.
326  """
327  if isinstance(index, numbers.Number):
328  index = self.getPatchIndexPairgetPatchIndexPair(index)
329  if (not 0 <= index[0] < self._numPatches_numPatches[0]) \
330  or (not 0 <= index[1] < self._numPatches_numPatches[1]):
331  raise IndexError("Patch index %s is not in range [0-%d, 0-%d]" %
332  (index, self._numPatches_numPatches[0]-1, self._numPatches_numPatches[1]-1))
333  innerMin = geom.Point2I(*[index[i] * self._patchInnerDimensions_patchInnerDimensions[i] for i in range(2)])
334  innerBBox = geom.Box2I(innerMin, self._patchInnerDimensions_patchInnerDimensions)
335  if not self._bbox.contains(innerBBox):
336  raise RuntimeError(
337  "Bug: patch index %s valid but inner bbox=%s not contained in tract bbox=%s" %
338  (index, innerBBox, self._bbox))
339  outerBBox = geom.Box2I(innerBBox)
340  outerBBox.grow(self.getPatchBordergetPatchBorder())
341  outerBBox.clip(self._bbox)
342  return PatchInfo(
343  index=index,
344  innerBBox=innerBBox,
345  outerBBox=outerBBox,
346  )
347 
349  """Get dimensions of inner region of the patches (all are the same)
350  """
351  return self._patchInnerDimensions_patchInnerDimensions
352 
353  def getTractOverlap(self):
354  """Get minimum overlap of adjacent sky tracts.
355  """
356  return self._tractOverlap_tractOverlap
357 
358  def getVertexList(self):
359  """Get list of ICRS sky coordinates of vertices that define the
360  boundary of the inner region.
361 
362  Notes
363  -----
364  **warning:** this is not a deep copy.
365  """
366  return self._vertexCoordList_vertexCoordList
367 
369  """Get inner on-sky region as a sphgeom.ConvexPolygon.
370  """
371  skyUnitVectors = [sp.getVector() for sp in self.getVertexListgetVertexList()]
372  return ConvexPolygon.convexHull(skyUnitVectors)
373 
375  """Get outer on-sky region as a sphgeom.ConvexPolygon
376  """
377  return makeSkyPolygonFromBBox(bbox=self.getBBoxgetBBox(), wcs=self.getWcsgetWcs())
378 
379  def getWcs(self):
380  """Get WCS of tract.
381 
382  Returns
383  -------
384  wcs : `lsst.afw.geom.SkyWcs`
385  The WCS of this tract
386  """
387  return self._wcs_wcs
388 
389  def __str__(self):
390  return "TractInfo(id=%s)" % (self._id_id,)
391 
392  def __repr__(self):
393  return "TractInfo(id=%s, ctrCoord=%s)" % (self._id_id, self._ctrCoord_ctrCoord.getVector())
394 
395  def __iter__(self):
396  xNum, yNum = self.getNumPatchesgetNumPatches()
397  for y in range(yNum):
398  for x in range(xNum):
399  yield self.getPatchInfogetPatchInfo((x, y))
400 
401  def __len__(self):
402  xNum, yNum = self.getNumPatchesgetNumPatches()
403  return xNum*yNum
404 
405  def __getitem__(self, index):
406  return self.getPatchInfogetPatchInfo(index)
407 
408  def contains(self, coord):
409  """Does this tract contain the coordinate?"""
410  try:
411  pixels = self.getWcsgetWcs().skyToPixel(coord)
413  # Point must be way off the tract
414  return False
415  return self.getBBoxgetBBox().contains(geom.Point2I(pixels))
416 
417 
419  """Information for a tract specified explicitly.
420 
421  A tract is placed at the explicitly defined coordinates, with the nominated
422  radius. The tracts are square (i.e., the radius is really a half-size).
423  """
424 
425  def __init__(self, ident, patchInnerDimensions, patchBorder, ctrCoord, radius, tractOverlap, wcs):
426  # We don't want TractInfo setting the bbox on the basis of vertices, but on the radius.
427  vertexList = []
428  self._radius_radius = radius
429  super(ExplicitTractInfo, self).__init__(ident, patchInnerDimensions, patchBorder, ctrCoord,
430  vertexList, tractOverlap, wcs)
431  # Shrink the box slightly to make sure the vertices are in the tract
432  bboxD = geom.BoxD(self.getBBoxgetBBox())
433  bboxD.grow(-0.001)
434  finalWcs = self.getWcsgetWcs()
435  self._vertexCoordList_vertexCoordList_vertexCoordList = finalWcs.pixelToSky(bboxD.getCorners())
436 
437  def _minimumBoundingBox(self, wcs):
438  """Calculate the minimum bounding box for the tract, given the WCS, and
439  the nominated radius.
440  """
441  bbox = geom.Box2D()
442  for i in range(4):
443  cornerCoord = self._ctrCoord_ctrCoord.offset(i*90*geom.degrees, self._radius_radius + self._tractOverlap_tractOverlap)
444  pixPos = wcs.skyToPixel(cornerCoord)
445  bbox.include(pixPos)
446  return bbox
A class representing an angle.
Definition: Angle.h:127
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
An integer coordinate rectangle.
Definition: Box.h:55
Reports arguments outside the domain of an operation.
Definition: Runtime.h:57
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
def __init__(self, ident, patchInnerDimensions, patchBorder, ctrCoord, radius, tractOverlap, wcs)
Definition: tractInfo.py:425
def _finalOrientation(self, bbox, wcs)
Definition: tractInfo.py:161
def getSequentialPatchIndex(self, patchInfo)
Definition: tractInfo.py:189
def findPatchList(self, coordList)
Definition: tractInfo.py:233
def __init__(self, id, patchInnerDimensions, patchBorder, ctrCoord, vertexCoordList, tractOverlap, wcs)
Definition: tractInfo.py:85
def getPatchInfo(self, index)
Definition: tractInfo.py:307
def __getitem__(self, index)
Definition: tractInfo.py:405
def getPatchIndexPair(self, sequentialIndex)
Definition: tractInfo.py:197
def _setupPatches(self, minBBox, wcs)
Definition: tractInfo.py:125
def findPatch(self, coord)
Definition: tractInfo.py:203
def _minimumBoundingBox(self, wcs)
Definition: tractInfo.py:101
def makeSkyPolygonFromBBox(bbox, wcs)
Definition: patchInfo.py:29