LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
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 numpy as np
26 
28 import lsst.geom as geom
29 from lsst.sphgeom import ConvexPolygon
30 
31 from .detail import makeSkyPolygonFromBBox, Index2D
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  tractBuilder : Subclass of `lsst.skymap.BaseTractBuilder`
42  Object used to compute patch geometry.
43  ctrCoord : `lsst.geom.SpherePoint`
44  ICRS sky coordinate of center of inner region of tract; also used as
45  the CRVAL for the WCS.
46  vertexCoordList : `list` of `lsst.geom.SpherePoint`
47  Vertices that define the boundaries of the inner region.
48  tractOverlap : `lsst.geom.Angle`
49  Minimum overlap between adjacent sky tracts; this defines the minimum
50  distance the tract extends beyond the inner region in all directions.
51  wcs : `lsst.afw.image.SkyWcs`
52  WCS for tract. The reference pixel will be shifted as required so that
53  the lower left-hand pixel (index 0,0) has pixel position 0.0, 0.0.
54 
55  Notes
56  -----
57  The tract is subdivided into rectangular patches. Each patch has the
58  following properties:
59 
60  - An inner region defined by an inner bounding box. The inner regions of
61  the patches exactly tile the tract, and all inner regions have the same
62  dimensions. The tract is made larger as required to make this work.
63 
64  - An outer region defined by an outer bounding box. The outer region
65  extends beyond the inner region by patchBorder pixels in all directions,
66  except there is no border at the edges of the tract.
67  Thus patches overlap each other but never extend off the tract.
68  If you do not want any overlap between adjacent patches then set
69  patchBorder to 0.
70 
71  - An index that consists of a pair of integers:
72 
73  * 0 <= x index < numPatches[0]
74 
75  * 0 <= y index < numPatches[1]
76 
77  Patch 0,0 is at the minimum corner of the tract bounding box.
78 
79  - It is not enforced that ctrCoord is the center of vertexCoordList, but
80  SkyMap relies on it.
81  """
82  def __init__(self, id, tractBuilder, ctrCoord, vertexCoordList, tractOverlap, wcs):
83  self._id_id = id
84  self._ctrCoord_ctrCoord = ctrCoord
85  self._vertexCoordList_vertexCoordList = tuple(vertexCoordList)
86  self._tractOverlap_tractOverlap = tractOverlap
87  self._tractBuilder_tractBuilder = tractBuilder
88 
89  minBBox = self._minimumBoundingBox_minimumBoundingBox(wcs)
90  initialBBox, self._numPatches_numPatches = self._tractBuilder_tractBuilder.setupPatches(minBBox, wcs)
91  self._bbox, self._wcs_wcs = self._finalOrientation_finalOrientation(initialBBox, wcs)
92 
93  def _minimumBoundingBox(self, wcs):
94  """Calculate the minimum bounding box for the tract, given the WCS.
95 
96  The bounding box is created in the frame of the supplied WCS,
97  so that it's OK if the coordinates are negative.
98 
99  We compute the bounding box that holds all the vertices and the
100  desired overlap.
101  """
102  minBBoxD = geom.Box2D()
103  halfOverlap = self._tractOverlap_tractOverlap / 2.0
104  for vertexCoord in self._vertexCoordList_vertexCoordList:
105  if self._tractOverlap_tractOverlap == 0:
106  minBBoxD.include(wcs.skyToPixel(vertexCoord))
107  else:
108  numAngles = 24
109  angleIncr = geom.Angle(360.0, geom.degrees) / float(numAngles)
110  for i in range(numAngles):
111  offAngle = angleIncr * i
112  offCoord = vertexCoord.offset(offAngle, halfOverlap)
113  pixPos = wcs.skyToPixel(offCoord)
114  minBBoxD.include(pixPos)
115  return minBBoxD
116 
117  def _finalOrientation(self, bbox, wcs):
118  """Determine the final orientation
119 
120  We offset everything so the lower-left corner is at 0,0
121  and compute the final Wcs.
122 
123  Parameters
124  ----------
125  bbox : `lsst.geom.Box2I`
126  Current bounding box.
127  wcs : `lsst.afw.geom.SkyWcs
128  Current Wcs.
129 
130  Returns
131  -------
132  finalBBox : `lsst.geom.Box2I`
133  Revised bounding box.
134  wcs : `lsst.afw.geom.SkyWcs`
135  Revised Wcs.
136  """
137  finalBBox = geom.Box2I(geom.Point2I(0, 0), bbox.getDimensions())
138  # shift the WCS by the same amount as the bbox; extra code is required
139  # because simply subtracting makes an Extent2I
140  pixPosOffset = geom.Extent2D(finalBBox.getMinX() - bbox.getMinX(),
141  finalBBox.getMinY() - bbox.getMinY())
142  wcs = wcs.copyAtShiftedPixelOrigin(pixPosOffset)
143  return finalBBox, wcs
144 
145  def getSequentialPatchIndex(self, patchInfo):
146  """Return a single integer that uniquely identifies
147  the given patch within this tract.
148 
149  Parameters
150  ----------
151  patchInfo : `lsst.skymap.PatchInfo`
152 
153  Returns
154  -------
155  sequentialIndex : `int`
156  """
157  return self._tractBuilder_tractBuilder.getSequentialPatchIndex(patchInfo)
158 
160  """Return a single integer that uniquely identifies
161  the patch index within the tract.
162 
163  Parameters
164  ----------
165  index : `lsst.skymap.Index2D`
166 
167  Returns
168  -------
169  sequentialIndex : `int`
170  """
171  return self._tractBuilder_tractBuilder.getSequentialPatchIndexFromPair(index)
172 
173  def getPatchIndexPair(self, sequentialIndex):
174  """Convert sequential index into patch index (x,y) pair.
175 
176  Parameters
177  ----------
178  sequentialIndex : `int`
179 
180  Returns
181  -------
182  x, y : `lsst.skymap.Index2D`
183  """
184  return self._tractBuilder_tractBuilder.getPatchIndexPair(sequentialIndex)
185 
186  def findPatch(self, coord):
187  """Find the patch containing the specified coord.
188 
189  Parameters
190  ----------
191  coord : `lsst.geom.SpherePoint`
192  ICRS sky coordinate to search for.
193 
194  Returns
195  -------
196  result : `lsst.skymap.PatchInfo`
197  PatchInfo of patch whose inner bbox contains the specified coord
198 
199  Raises
200  ------
201  LookupError
202  If coord is not in tract or we cannot determine the
203  pixel coordinate (which likely means the coord is off the tract).
204  """
205  try:
206  pixel = self.wcswcs.skyToPixel(coord)
208  # Point must be way off the tract
209  raise LookupError("Unable to determine pixel position for coordinate %s" % (coord,))
210  pixelInd = geom.Point2I(pixel)
211  if not self._bbox.contains(pixelInd):
212  raise LookupError("coord %s is not in tract %s" % (coord, self.tract_idtract_id))
213  # This should probably be the same as above because we only
214  # care about the INNER dimensions.
215  patchInd = tuple(int(pixelInd[i]/self.patch_inner_dimensionspatch_inner_dimensions[i]) for i in range(2))
216  return self.getPatchInfogetPatchInfo(patchInd)
217 
218  def findPatchList(self, coordList):
219  """Find patches containing the specified list of coords.
220 
221  Parameters
222  ----------
223  coordList : `list` of `lsst.geom.SpherePoint`
224  ICRS sky coordinates to search for.
225 
226  Returns
227  -------
228  result : `list` of `lsst.skymap.PatchInfo`
229  List of PatchInfo for patches that contain, or may contain, the
230  specified region. The list will be empty if there is no overlap.
231 
232  Notes
233  -----
234  **Warning:**
235 
236  - This may give incorrect answers on regions that are larger than a
237  tract.
238 
239  - This uses a naive algorithm that may find some patches that do not
240  overlap the region (especially if the region is not a rectangle
241  aligned along patch x,y).
242  """
243  box2D = geom.Box2D()
244  for coord in coordList:
245  try:
246  pixelPos = self.wcswcs.skyToPixel(coord)
248  # the point is so far off the tract that its pixel position cannot be computed
249  continue
250  box2D.include(pixelPos)
251  bbox = geom.Box2I(box2D)
252  bbox.grow(self.getPatchBordergetPatchBorder())
253  bbox.clip(self._bbox)
254  if bbox.isEmpty():
255  return ()
256 
257  llPatchInd = tuple(int(bbox.getMin()[i]/self.patch_inner_dimensionspatch_inner_dimensions[i]) for i in range(2))
258  urPatchInd = tuple(int(bbox.getMax()[i]/self.patch_inner_dimensionspatch_inner_dimensions[i]) for i in range(2))
259  return tuple(self.getPatchInfogetPatchInfo((xInd, yInd))
260  for xInd in range(llPatchInd[0], urPatchInd[0]+1)
261  for yInd in range(llPatchInd[1], urPatchInd[1]+1))
262 
263  def getBBox(self):
264  """Get bounding box of tract (as an geom.Box2I)
265  """
266  return geom.Box2I(self._bbox)
267 
268  bbox = property(getBBox)
269 
270  def getCtrCoord(self):
271  """Get ICRS sky coordinate of center of tract
272  (as an lsst.geom.SpherePoint)
273  """
274  return self._ctrCoord_ctrCoord
275 
276  ctr_coord = property(getCtrCoord)
277 
278  def getId(self):
279  """Get ID of tract
280  """
281  return self._id_id
282 
283  tract_id = property(getId)
284 
285  def getNumPatches(self):
286  """Get the number of patches in x, y.
287 
288  Returns
289  -------
290  result : `lsst.skymap.Index2D`
291  The number of patches in x, y
292  """
293  return self._numPatches_numPatches
294 
295  num_patches = property(getNumPatches)
296 
297  def getPatchBorder(self):
298  return self._tractBuilder_tractBuilder.getPatchBorder()
299 
300  patch_border = property(getPatchBorder)
301 
302  def getPatchInfo(self, index):
303  """Return information for the specified patch.
304 
305  Parameters
306  ----------
307  index : `typing.NamedTuple` ['x': `int`, 'y': `int`]
308  Index of patch, as a pair of ints;
309  or a sequential index as returned by getSequentialPatchIndex;
310  negative values are not supported.
311 
312  Returns
313  -------
314  result : `lsst.skymap.PatchInfo`
315  The patch info for that index.
316 
317  Raises
318  ------
319  IndexError
320  If index is out of range.
321  """
322  return self._tractBuilder_tractBuilder.getPatchInfo(index, self._wcs_wcs)
323 
325  """Get dimensions of inner region of the patches (all are the same)
326  """
327  return self._tractBuilder_tractBuilder.getPatchInnerDimensions()
328 
329  patch_inner_dimensions = property(getPatchInnerDimensions)
330 
331  def getTractOverlap(self):
332  """Get minimum overlap of adjacent sky tracts.
333  """
334  return self._tractOverlap_tractOverlap
335 
336  tract_overlap = property(getTractOverlap)
337 
338  def getVertexList(self):
339  """Get list of ICRS sky coordinates of vertices that define the
340  boundary of the inner region.
341 
342  Notes
343  -----
344  **warning:** this is not a deep copy.
345  """
346  return self._vertexCoordList_vertexCoordList
347 
348  vertex_list = property(getVertexList)
349 
351  """Get inner on-sky region as a sphgeom.ConvexPolygon.
352  """
353  skyUnitVectors = [sp.getVector() for sp in self.getVertexListgetVertexList()]
354  return ConvexPolygon.convexHull(skyUnitVectors)
355 
356  inner_sky_polygon = property(getInnerSkyPolygon)
357 
359  """Get outer on-sky region as a sphgeom.ConvexPolygon
360  """
361  return makeSkyPolygonFromBBox(bbox=self.getBBoxgetBBox(), wcs=self.getWcsgetWcs())
362 
363  outer_sky_polygon = property(getOuterSkyPolygon)
364 
365  def getWcs(self):
366  """Get WCS of tract.
367 
368  Returns
369  -------
370  wcs : `lsst.afw.geom.SkyWcs`
371  The WCS of this tract
372  """
373  return self._wcs_wcs
374 
375  wcs = property(getWcs)
376 
377  def __str__(self):
378  return "TractInfo(id=%s)" % (self._id_id,)
379 
380  def __repr__(self):
381  return "TractInfo(id=%s, ctrCoord=%s)" % (self._id_id, self._ctrCoord_ctrCoord.getVector())
382 
383  def __iter__(self):
384  xNum, yNum = self.getNumPatchesgetNumPatches()
385  for y in range(yNum):
386  for x in range(xNum):
387  yield self.getPatchInfogetPatchInfo(Index2D(x=x, y=y))
388 
389  def __len__(self):
390  xNum, yNum = self.getNumPatchesgetNumPatches()
391  return xNum*yNum
392 
393  def __getitem__(self, index):
394  return self.getPatchInfogetPatchInfo(index)
395 
396  def contains(self, coord):
397  """Does this tract contain the coordinate?"""
398  try:
399  pixel = self.getWcsgetWcs().skyToPixel(coord)
401  # Point must be way off the tract
402  return False
403  if not np.isfinite(pixel.getX()) or not np.isfinite(pixel.getY()):
404  # Point is definitely off the tract
405  return False
406  return self.getBBoxgetBBox().contains(geom.Point2I(pixel))
407 
408 
410  """Information for a tract specified explicitly.
411 
412  A tract is placed at the explicitly defined coordinates, with the nominated
413  radius. The tracts are square (i.e., the radius is really a half-size).
414  """
415  def __init__(self, id, tractBuilder, ctrCoord, radius, tractOverlap, wcs):
416  # We don't want TractInfo setting the bbox on the basis of vertices, but on the radius.
417  vertexList = []
418  self._radius_radius = radius
419  super(ExplicitTractInfo, self).__init__(id, tractBuilder, ctrCoord,
420  vertexList, tractOverlap, wcs)
421  # Shrink the box slightly to make sure the vertices are in the tract
422  bboxD = geom.BoxD(self.getBBoxgetBBox())
423  bboxD.grow(-0.001)
424  finalWcs = self.getWcsgetWcs()
425  self._vertexCoordList_vertexCoordList_vertexCoordList = finalWcs.pixelToSky(bboxD.getCorners())
426 
427  def _minimumBoundingBox(self, wcs):
428  """Calculate the minimum bounding box for the tract, given the WCS, and
429  the nominated radius.
430  """
431  bbox = geom.Box2D()
432  for i in range(4):
433  cornerCoord = self._ctrCoord_ctrCoord.offset(i*90*geom.degrees, self._radius_radius + self._tractOverlap_tractOverlap)
434  pixPos = wcs.skyToPixel(cornerCoord)
435  bbox.include(pixPos)
436  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, id, tractBuilder, ctrCoord, radius, tractOverlap, wcs)
Definition: tractInfo.py:415
def _finalOrientation(self, bbox, wcs)
Definition: tractInfo.py:117
def getSequentialPatchIndex(self, patchInfo)
Definition: tractInfo.py:145
def findPatchList(self, coordList)
Definition: tractInfo.py:218
def getPatchInfo(self, index)
Definition: tractInfo.py:302
def __getitem__(self, index)
Definition: tractInfo.py:393
def __init__(self, id, tractBuilder, ctrCoord, vertexCoordList, tractOverlap, wcs)
Definition: tractInfo.py:82
def getSequentialPatchIndexFromPair(self, index)
Definition: tractInfo.py:159
def getPatchIndexPair(self, sequentialIndex)
Definition: tractInfo.py:173
def findPatch(self, coord)
Definition: tractInfo.py:186
def _minimumBoundingBox(self, wcs)
Definition: tractInfo.py:93
def makeSkyPolygonFromBBox(bbox, wcs)
Definition: utils.py:61