LSSTApplications  11.0-24-g0a022a1,14.0+77,15.0,15.0+1
LSSTDataManagementBasePackage
quadsphere.py
Go to the documentation of this file.
1 from __future__ import division
2 from builtins import map
3 from builtins import range
4 from builtins import object
5 #
6 # LSST Data Management System
7 # Copyright 2008, 2009, 2010 LSST Corporation.
8 #
9 # This product includes software developed by the
10 # LSST Project (http://www.lsst.org/).
11 #
12 # This program is free software: you can redistribute it and/or modify
13 # it under the terms of the GNU General Public License as published by
14 # the Free Software Foundation, either version 3 of the License, or
15 # (at your option) any later version.
16 #
17 # This program is distributed in the hope that it will be useful,
18 # but WITHOUT ANY WARRANTY; without even the implied warranty of
19 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 # GNU General Public License for more details.
21 #
22 # You should have received a copy of the LSST License Statement and
23 # the GNU General Public License along with this program. If not,
24 # see <http://www.lsstcorp.org/LegalNotices/>.
25 #
26 
27 import math
28 import numbers
29 
30 import lsst.pex.policy as pexPolicy
31 import lsst.geom as geom
32 
33 
34 # Methods for rotating vectors around the x, y, z and -x, -y, -z axes
35 def _rotX(v, sin, cos):
36  return (v[0], v[1] * cos - v[2] * sin, v[1] * sin + v[2] * cos)
37 
38 
39 def _rotY(v, sin, cos):
40  return (v[0] * cos + v[2] * sin, v[1], -v[0] * sin + v[2] * cos)
41 
42 
43 def _rotZ(v, sin, cos):
44  return (v[0] * cos - v[1] * sin, v[0] * sin + v[1] * cos, v[2])
45 
46 
47 def _rotNX(v, sin, cos):
48  return _rotX(v, -sin, cos)
49 
50 
51 def _rotNY(v, sin, cos):
52  return _rotY(v, -sin, cos)
53 
54 
55 def _rotNZ(v, sin, cos):
56  return _rotZ(v, -sin, cos)
57 
58 
59 class QuadSpherePixelization(object):
60  """A quad-sphere based pixelization of the unit sphere. Each of the 6
61  root pixels are divided into an R by R grid of sky-pixels, where R
62  is the resolution (a positive integer). Sky-pixel identifiers are
63  contiguous integers in range [0, 6 * R**2), and are ordered such that
64  enumerating the pixels 0, 1, ..., 6 * R**2 results in a roughly spiral
65  traversal (from south to north pole) of the unit sphere. In particular,
66  iterating over pixels in this order guarantees that when arriving at
67  pixel P:
68 
69  - all pixels with id P' < P have been visited
70  - all neighbors of pixels with id P' < P - 8R have been visited
71 
72  Each pixel may optionally be expanded by a padding angle A, such that
73  points on the unit sphere outside of a pixel but within angular separation
74  A of the fiducial pixel boundaries are also considered to be part of that
75  pixel. Pixel (and padded pixel) edges are great circles on the unit
76  sphere.
77 
78  This class provides methods for computing fiducial or padded geometric
79  representations of a sky-pixel - spherical convex polygons in both
80  cases. It also supports finding the list of sky-pixels intersecting
81  an input spherical convex polygon, as well as determining the neighbors
82  and center of a sky-pixel.
83 
84  TODO: Right now, pixels are defined by a simple tangent plane projection
85  of each cube face onto the sphere, but pixel boundaries could be adjusted
86  to produce equal area pixels.
87  """
88 
89  def __init__(self, resolution, paddingRad):
90  """Creates a new quad-sphere sky pixelisation.
91 
92  resolution: the number of pixels along the x and y axes of each root
93  pixel.
94 
95  paddingRad: the angular separation (rad) by which fiducial sky-pixels
96  are to be padded.
97  """
98  if not isinstance(resolution, numbers.Integral):
99  raise TypeError(
100  'Quad-sphere resolution parameter must be an integer')
101  if resolution < 3 or resolution > 16384:
102  raise RuntimeError(
103  'Quad-sphere resolution must be in range [3, 16384]')
104  if not isinstance(paddingRad, float):
105  raise TypeError(
106  'Quad-sphere pixel padding angle must be a float')
107  if paddingRad < 0.0 or paddingRad >= math.pi * 0.25:
108  raise RuntimeError(
109  'Quad-sphere pixel padding radius must be in range [0, 45) deg')
110  R = resolution
111  self.resolution = R
112  self.padding = paddingRad
113  x, y, z = (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)
114  nx, ny, nz = (-1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, -1.0)
115  # center of each root pixel
116  self.center = [z, x, y, nx, ny, nz]
117  # x basis vector of each root pixel
118  self.x = [ny, y, nx, ny, x, ny]
119  # y basis vector of each root pixel
120  self.y = [x, z, z, z, z, nx]
121  # functions for rotating vectors in direction of increasing x (per root pixel)
122  self.xrot = [_rotX, _rotZ, _rotZ, _rotZ, _rotZ, _rotNX]
123  # functions for rotating vectors in direction of increasing y (per root pixel)
124  self.yrot = [_rotY, _rotNY, _rotX, _rotY, _rotNX, _rotY]
125  # Compute fiducial and splitting x/y planes for each root pixel
126  self.xplane = []
127  self.yplane = []
128  sp = math.sin(self.padding)
129  for root in range(6):
130  xplanes = []
131  yplanes = []
132  c, x, y = self.center[root], self.x[root], self.y[root]
133  for i in range(R + 1):
134  xfp = self._fiducialXPlane(root, i)
135  yfp = self._fiducialYPlane(root, i)
136  f = 2.0 * float(i) / float(R) - 1.0
137  thetaX = math.radians(0.5 * geom.cartesianAngularSep(
138  (c[0] + f * x[0] + y[0],
139  c[1] + f * x[1] + y[1],
140  c[2] + f * x[2] + y[2]),
141  (c[0] + f * x[0] - y[0],
142  c[1] + f * x[1] - y[1],
143  c[2] + f * x[2] - y[2])))
144  thetaY = math.radians(0.5 * geom.cartesianAngularSep(
145  (c[0] + x[0] + f * y[0],
146  c[1] + x[1] + f * y[1],
147  c[2] + x[2] + f * y[2]),
148  (c[0] - x[0] + f * y[0],
149  c[1] - x[1] + f * y[1],
150  c[2] - x[2] + f * y[2])))
151  sinrx = sp / math.cos(thetaX)
152  sinry = sp / math.cos(thetaY)
153  cosrx = math.sqrt(1.0 - sinrx * sinrx)
154  cosry = math.sqrt(1.0 - sinry * sinry)
155  if i == 0:
156  xlp = None
157  ybp = None
158  else:
159  xlp = self.xrot[root](xfp, sinrx, cosrx)
160  ybp = self.yrot[root](yfp, sinry, cosry)
161  xlp = (-xlp[0], -xlp[1], -xlp[2])
162  ybp = (-ybp[0], -ybp[1], -ybp[2])
163  if i == R:
164  xrp = None
165  ytp = None
166  else:
167  xrp = self.xrot[root](xfp, -sinrx, cosrx)
168  ytp = self.yrot[root](yfp, -sinry, cosry)
169  xplanes.append((xfp, xlp, xrp))
170  yplanes.append((yfp, ybp, ytp))
171  self.xplane.append(xplanes)
172  self.yplane.append(yplanes)
173  # Corner pixel neighbors.
174  self.cornerNeighbors = (
175  # root pixel 0
176  (
177  # x, y = 0, 0
178  (self.id(0, 1, 0),
179  self.id(0, 1, 1),
180  self.id(0, 0, 1),
181  self.id(2, R - 1, R - 1),
182  self.id(2, R - 2, R - 1),
183  self.id(3, 0, R - 1),
184  self.id(3, 1, R - 1)),
185  # x, y = 0, R-1
186  (self.id(0, 1, R - 1),
187  self.id(0, 1, R - 2),
188  self.id(0, 0, R - 2),
189  self.id(1, R - 1, R - 1),
190  self.id(1, R - 2, R - 1),
191  self.id(2, 0, R - 1),
192  self.id(2, 1, R - 1)),
193  # x, y = R-1, 0
194  (self.id(0, R - 2, 0),
195  self.id(0, R - 2, 1),
196  self.id(0, R - 1, 1),
197  self.id(3, R - 2, R - 1),
198  self.id(3, R - 1, R - 1),
199  self.id(4, 0, R - 1),
200  self.id(4, 1, R - 1)),
201  # x, y = R-1, R-1
202  (self.id(0, R - 2, R - 1),
203  self.id(0, R - 2, R - 2),
204  self.id(0, R - 1, R - 2),
205  self.id(1, 0, R - 1),
206  self.id(1, 1, R - 1),
207  self.id(4, R - 2, R - 1),
208  self.id(4, R - 1, R - 1))
209  ),
210  # root pixel 1
211  (
212  # x, y = 0, 0
213  (self.id(1, 1, 0),
214  self.id(1, 1, 1),
215  self.id(1, 0, 1),
216  self.id(4, R - 1, 0),
217  self.id(4, R - 1, 1),
218  self.id(5, R - 2, 0),
219  self.id(5, R - 1, 0)),
220  # x, y = 0, R-1
221  (self.id(1, 1, R - 1),
222  self.id(1, 1, R - 2),
223  self.id(1, 0, R - 2),
224  self.id(4, R - 1, R - 2),
225  self.id(4, R - 1, R - 1),
226  self.id(0, R - 2, R - 1),
227  self.id(0, R - 1, R - 1)),
228  # x, y = R-1, 0
229  (self.id(1, R - 2, 0),
230  self.id(1, R - 2, 1),
231  self.id(1, R - 1, 1),
232  self.id(2, 0, 0),
233  self.id(2, 0, 1),
234  self.id(5, 0, 0),
235  self.id(5, 1, 0)),
236  # x, y = R-1, R-1
237  (self.id(1, R - 2, R - 1),
238  self.id(1, R - 2, R - 2),
239  self.id(1, R - 1, R - 2),
240  self.id(2, 0, R - 2),
241  self.id(2, 0, R - 1),
242  self.id(0, 0, R - 1),
243  self.id(0, 1, R - 1))
244  ),
245  # root pixel 2
246  (
247  # x, y = 0, 0
248  (self.id(2, 1, 0),
249  self.id(2, 1, 1),
250  self.id(2, 0, 1),
251  self.id(1, R - 1, 0),
252  self.id(1, R - 1, 1),
253  self.id(5, 0, 0),
254  self.id(5, 0, 1)),
255  # x, y = 0, R-1
256  (self.id(2, 1, R - 1),
257  self.id(2, 1, R - 2),
258  self.id(2, 0, R - 2),
259  self.id(1, R - 1, R - 2),
260  self.id(1, R - 1, R - 1),
261  self.id(0, 0, R - 2),
262  self.id(0, 0, R - 1)),
263  # x, y = R-1, 0
264  (self.id(2, R - 2, 0),
265  self.id(2, R - 2, 1),
266  self.id(2, R - 1, 1),
267  self.id(3, 0, 0),
268  self.id(3, 0, 1),
269  self.id(5, 0, R - 2),
270  self.id(5, 0, R - 1)),
271  # x, y = R-1, R-1
272  (self.id(2, R - 2, R - 1),
273  self.id(2, R - 2, R - 2),
274  self.id(2, R - 1, R - 2),
275  self.id(3, 0, R - 2),
276  self.id(3, 0, R - 1),
277  self.id(0, 0, 0),
278  self.id(0, 0, 1))
279  ),
280  # root pixel 3
281  (
282  # x, y = 0, 0
283  (self.id(3, 1, 0),
284  self.id(3, 1, 1),
285  self.id(3, 0, 1),
286  self.id(2, R - 1, 0),
287  self.id(2, R - 1, 1),
288  self.id(5, 0, R - 1),
289  self.id(5, 1, R - 1)),
290  # x, y = 0, R-1
291  (self.id(3, 1, R - 1),
292  self.id(3, 1, R - 2),
293  self.id(3, 0, R - 2),
294  self.id(2, R - 1, R - 2),
295  self.id(2, R - 1, R - 1),
296  self.id(0, 0, 0),
297  self.id(0, 1, 0)),
298  # x, y = R-1, 0
299  (self.id(3, R - 2, 0),
300  self.id(3, R - 2, 1),
301  self.id(3, R - 1, 1),
302  self.id(4, 0, 0),
303  self.id(4, 0, 1),
304  self.id(5, R - 2, R - 1),
305  self.id(5, R - 1, R - 1)),
306  # x, y = R-1, R-1
307  (self.id(3, R - 2, R - 1),
308  self.id(3, R - 2, R - 2),
309  self.id(3, R - 1, R - 2),
310  self.id(4, 0, R - 2),
311  self.id(4, 0, R - 1),
312  self.id(0, R - 2, 0),
313  self.id(0, R - 1, 0))
314  ),
315  # root pixel 4
316  (
317  # x, y = 0, 0
318  (self.id(4, 1, 0),
319  self.id(4, 1, 1),
320  self.id(4, 0, 1),
321  self.id(3, R - 1, 0),
322  self.id(3, R - 1, 1),
323  self.id(5, R - 1, R - 2),
324  self.id(5, R - 1, R - 1)),
325  # x, y = 0, R-1
326  (self.id(4, 1, R - 1),
327  self.id(4, 1, R - 2),
328  self.id(4, 0, R - 2),
329  self.id(3, R - 1, R - 2),
330  self.id(3, R - 1, R - 1),
331  self.id(0, R - 1, 0),
332  self.id(0, R - 1, 1)),
333  # x, y = R-1, 0
334  (self.id(4, R - 2, 0),
335  self.id(4, R - 2, 1),
336  self.id(4, R - 1, 1),
337  self.id(1, 0, 0),
338  self.id(1, 0, 1),
339  self.id(5, R - 1, 0),
340  self.id(5, R - 1, 1)),
341  # x, y = R-1, R-1
342  (self.id(4, R - 2, R - 1),
343  self.id(4, R - 2, R - 2),
344  self.id(4, R - 1, R - 2),
345  self.id(1, 0, R - 2),
346  self.id(1, 0, R - 1),
347  self.id(0, R - 1, R - 2),
348  self.id(0, R - 1, R - 1))
349  ),
350  # root pixel 5
351  (
352  # x, y = 0, 0
353  (self.id(5, 1, 0),
354  self.id(5, 1, 1),
355  self.id(5, 0, 1),
356  self.id(1, R - 2, 0),
357  self.id(1, R - 1, 0),
358  self.id(2, 0, 0),
359  self.id(2, 1, 0)),
360  # x, y = 0, R-1
361  (self.id(5, 1, R - 1),
362  self.id(5, 1, R - 2),
363  self.id(5, 0, R - 2),
364  self.id(2, R - 2, 0),
365  self.id(2, R - 1, 0),
366  self.id(3, 0, 0),
367  self.id(3, 1, 0)),
368  # x, y = R-1, 0
369  (self.id(5, R - 2, 0),
370  self.id(5, R - 2, 1),
371  self.id(5, R - 1, 1),
372  self.id(1, 0, 0),
373  self.id(1, 1, 0),
374  self.id(4, R - 2, 0),
375  self.id(4, R - 1, 0)),
376  # x, y = R-1, R-1
377  (self.id(5, R - 2, R - 1),
378  self.id(5, R - 2, R - 2),
379  self.id(5, R - 1, R - 2),
380  self.id(3, R - 2, 0),
381  self.id(3, R - 1, 0),
382  self.id(4, 0, 0),
383  self.id(4, 1, 0))
384  ),
385  )
386 
387  def getGeometry(self, pixelId, fiducial=False):
388  """Returns a spherical convex polygon corresponding to the fiducial
389  or padded boundaries of the sky-pixel with the specified id.
390  """
391  root, ix, iy = self.coords(pixelId)
392  xl = 2.0 * float(ix) / float(self.resolution) - 1.0
393  xr = 2.0 * float(ix + 1) / float(self.resolution) - 1.0
394  yb = 2.0 * float(iy) / float(self.resolution) - 1.0
395  yt = 2.0 * float(iy + 1) / float(self.resolution) - 1.0
396  c, x, y = self.center[root], self.x[root], self.y[root]
397  v = list(map(geom.normalize, [(c[0] + xl * x[0] + yb * y[0],
398  c[1] + xl * x[1] + yb * y[1],
399  c[2] + xl * x[2] + yb * y[2]),
400  (c[0] + xr * x[0] + yb * y[0],
401  c[1] + xr * x[1] + yb * y[1],
402  c[2] + xr * x[2] + yb * y[2]),
403  (c[0] + xr * x[0] + yt * y[0],
404  c[1] + xr * x[1] + yt * y[1],
405  c[2] + xr * x[2] + yt * y[2]),
406  (c[0] + xl * x[0] + yt * y[0],
407  c[1] + xl * x[1] + yt * y[1],
408  c[2] + xl * x[2] + yt * y[2])]))
409  if not fiducial and self.padding > 0.0:
410  # Determine angles by which edge planes must be rotated outwards
411  sp = math.sin(self.padding)
412  theta = [0.5 * geom.cartesianAngularSep(x[0], x[1]) for
413  x in ((v[0], v[3]), (v[1], v[0]), (v[2], v[1]), (v[3], v[2]))]
414  sina = [sp / math.cos(math.radians(x)) for x in theta]
415  cosa = [math.sqrt(1.0 - x * x) for x in sina]
416  # find plane equations of fiducial pixel boundaries
417  xlp = self.xplane[root][ix][0]
418  ybp = self.yplane[root][iy][0]
419  xrp = self.xplane[root][ix + 1][0]
420  ytp = self.yplane[root][iy + 1][0]
421  # rotate edge planes outwards
422  xlp = self.xrot[root](xlp, -sina[0], cosa[0])
423  ybp = self.yrot[root](ybp, -sina[1], cosa[1])
424  xrp = self.xrot[root](xrp, sina[2], cosa[2])
425  ytp = self.yrot[root](ytp, sina[3], cosa[3])
426  # intersect rotated planes to find vertices of padded sky-pixel
427  v = list(map(geom.normalize, [geom.cross(xlp, ybp),
428  geom.cross(xrp, ybp),
429  geom.cross(xrp, ytp),
430  geom.cross(xlp, ytp)]))
431  return geom.SphericalConvexPolygon(v)
432 
433  def getCenter(self, pixelId):
434  """Returns the center of a sky-pixel as a unit cartesian 3-vector.
435  """
436  root, ix, iy = self.coords(pixelId)
437  xc = 2.0 * (float(ix) + 0.5) / float(self.resolution) - 1.0
438  yc = 2.0 * (float(iy) + 0.5) / float(self.resolution) - 1.0
439  c, x, y = self.center[root], self.x[root], self.y[root]
440  return geom.normalize((c[0] + xc * x[0] + yc * y[0],
441  c[1] + xc * x[1] + yc * y[1],
442  c[2] + xc * x[2] + yc * y[2]))
443 
444  def getNeighbors(self, pixelId):
445  """Returns a list of ids for sky-pixels adjacent to specified pixel.
446  """
447  root, ix, iy = self.coords(pixelId)
448  R = self.resolution
449  if ix == 0:
450  # left edge
451  if iy == 0:
452  return self.cornerNeighbors[root][0]
453  elif iy == R - 1:
454  return self.cornerNeighbors[root][1]
455  elif root == 0:
456  neighbors = (self.id(2, R - iy - 2, R - 1),
457  self.id(2, R - iy - 1, R - 1),
458  self.id(2, R - iy, R - 1))
459  elif root == 5:
460  neighbors = (self.id(2, iy - 1, 0),
461  self.id(2, iy, 0),
462  self.id(2, iy + 1, 0))
463  else:
464  r = root - 1
465  if r == 0:
466  r = 4
467  neighbors = (self.id(r, R - 1, iy - 1),
468  self.id(r, R - 1, iy),
469  self.id(r, R - 1, iy + 1))
470  return neighbors + (self.id(root, 1, iy - 1),
471  self.id(root, 1, iy),
472  self.id(root, 1, iy + 1),
473  self.id(root, 0, iy - 1),
474  self.id(root, 0, iy + 1))
475  elif ix == R - 1:
476  # right edge
477  if iy == 0:
478  return self.cornerNeighbors[root][2]
479  elif iy == R - 1:
480  return self.cornerNeighbors[root][3]
481  elif root == 0:
482  neighbors = (self.id(4, iy - 1, R - 1),
483  self.id(4, iy, R - 1),
484  self.id(4, iy + 1, R - 1))
485  elif root == 5:
486  neighbors = (self.id(4, R - iy - 2, 0),
487  self.id(4, R - iy - 1, 0),
488  self.id(4, R - iy, 0))
489  else:
490  r = root + 1
491  if r == 5:
492  r = 1
493  neighbors = (self.id(r, 0, iy - 1),
494  self.id(r, 0, iy),
495  self.id(r, 0, iy + 1))
496  return neighbors + (self.id(root, R - 2, iy - 1),
497  self.id(root, R - 2, iy),
498  self.id(root, R - 2, iy + 1),
499  self.id(root, R - 1, iy - 1),
500  self.id(root, R - 1, iy + 1))
501  elif iy == 0:
502  # bottom edge
503  if root == 0:
504  neighbors = (self.id(3, ix - 1, R - 1),
505  self.id(3, ix, R - 1),
506  self.id(3, ix + 1, R - 1))
507  elif root == 1:
508  neighbors = (self.id(5, R - ix - 2, 0),
509  self.id(5, R - ix - 1, 0),
510  self.id(5, R - ix, 0))
511  elif root == 2:
512  neighbors = (self.id(5, 0, ix - 1),
513  self.id(5, 0, ix),
514  self.id(5, 0, ix + 1))
515  elif root == 3:
516  neighbors = (self.id(5, ix - 1, R - 1),
517  self.id(5, ix, R - 1),
518  self.id(5, ix + 1, R - 1))
519  elif root == 4:
520  neighbors = (self.id(5, R - 1, R - ix - 2),
521  self.id(5, R - 1, R - ix - 1),
522  self.id(5, R - 1, R - ix))
523  else:
524  neighbors = (self.id(1, R - ix - 2, 0),
525  self.id(1, R - ix - 1, 0),
526  self.id(1, R - ix, 0))
527  return neighbors + (self.id(root, ix - 1, 1),
528  self.id(root, ix, 1),
529  self.id(root, ix + 1, 1),
530  self.id(root, ix - 1, 0),
531  self.id(root, ix + 1, 0))
532  elif iy == R - 1:
533  # top edge
534  if root == 0:
535  neighbors = (self.id(1, R - ix - 2, R - 1),
536  self.id(1, R - ix - 1, R - 1),
537  self.id(1, R - ix, R - 1))
538  elif root == 1:
539  neighbors = (self.id(0, R - ix - 2, R - 1),
540  self.id(0, R - ix - 1, R - 1),
541  self.id(0, R - ix, R - 1))
542  elif root == 2:
543  neighbors = (self.id(0, 0, R - ix - 2),
544  self.id(0, 0, R - ix - 1),
545  self.id(0, 0, R - ix))
546  elif root == 3:
547  neighbors = (self.id(0, ix - 1, 0),
548  self.id(0, ix, 0),
549  self.id(0, ix + 1, 0))
550  elif root == 4:
551  neighbors = (self.id(0, R - 1, ix - 1),
552  self.id(0, R - 1, ix),
553  self.id(0, R - 1, ix + 1))
554  else:
555  neighbors = (self.id(3, ix - 1, 0),
556  self.id(3, ix, 0),
557  self.id(3, ix + 1, 0))
558  return neighbors + (self.id(root, ix - 1, R - 2),
559  self.id(root, ix, R - 2),
560  self.id(root, ix + 1, R - 2),
561  self.id(root, ix - 1, R - 1),
562  self.id(root, ix + 1, R - 1))
563  # interior pixel
564  return (self.id(root, ix - 1, iy - 1),
565  self.id(root, ix, iy - 1),
566  self.id(root, ix + 1, iy - 1),
567  self.id(root, ix - 1, iy),
568  self.id(root, ix + 1, iy),
569  self.id(root, ix - 1, iy + 1),
570  self.id(root, ix, iy + 1),
571  self.id(root, ix + 1, iy + 1))
572 
573  def intersect(self, polygon):
574  """Returns a list of ids for sky-pixels intersecting the given
575  polygon. Intersection is computed against padded sky-pixels.
576  """
577  if not isinstance(polygon, geom.SphericalConvexPolygon):
578  raise TypeError('Expecting a SphericalConvexPolygon as input ' +
579  'to the sky-pixel intersection computation')
580  pixels = []
581  R = self.resolution
582  for root in range(6):
583  # clip against padded root pixel edge planes
584  poly = polygon
585  for p in (self.xplane[root][0][2], self.xplane[root][R][1],
586  self.yplane[root][0][2], self.yplane[root][R][1]):
587  poly = poly.clip(p)
588  if poly is None:
589  break
590  if poly is None:
591  # polygon doesn't intersect this root triangle
592  continue
593  self._intersect(pixels, poly, root, (0, R, 0, R))
594  return pixels
595 
596  def __len__(self):
597  """Returns the total number of sky-pixels in this pixelization.
598  """
599  return 6 * self.resolution * self.resolution
600 
601  def __iter__(self):
602  """Returns an iterator over the sky-pixel ids of all the pixels
603  in this pixelization.
604  """
605  return iter(range(len(self)))
606 
607  def __repr__(self):
608  return ''.join([self.__class__.__name__, '(',
609  repr(self.resolution), ', ', repr(self.padding), ')'])
610 
611  def pixel(self, theta, phi):
612  """Maps from spherical coordinates (longitude and latitude
613  angles theta and phi, both in radians) to a sky-pixel id.
614  """
615  R = self.resolution
616  root = int(math.fmod(0.5 + 2.0 * theta / math.pi, 4.0))
617  theta1 = theta - 0.5 * math.pi * root
618  root += 1
619  tanPhi = math.tan(phi)
620  y = tanPhi / math.cos(theta1)
621  if y > 1:
622  root = 0
623  x = -math.sin(theta) / tanPhi
624  y = math.cos(theta) / tanPhi
625  elif y < -1:
626  root = 5
627  x = math.sin(theta) / tanPhi
628  y = math.cos(theta) / tanPhi
629  else:
630  x = math.tan(theta1)
631  x = int(R * 0.5 * (x + 1.0))
632  y = int(R * 0.5 * (y + 1.0))
633  if x >= R:
634  x = R - 1
635  if y >= R:
636  y = R - 1
637  return self.id(root, x, y)
638 
639  def id(self, root, x, y):
640  """Maps from a root pixel number and x, y pixel coordinates
641  to a sky-pixel id.
642  """
643  if not all(isinstance(p, numbers.Integral) for p in (root, x, y)):
644  raise TypeError('Pixel coordinates must all be an integer')
645  if root < 0 or root > 5:
646  raise RuntimeError('root pixel number must be in range [0,6)')
647  R = self.resolution
648  if x < 0 or x >= R or y < 0 or y >= R:
649  raise RuntimeError('x and y pixel coordinates must ' +
650  'be in range [0,%d)' % R)
651  if root > 0 and root < 5:
652  # easy equatorial pixel case
653  return R ** 2 + 4 * R * y + (root - 1) * R + x
654  # assign ids to pixels in expanding concentric squares
655  dx = x - 0.5 * (R - 1)
656  dy = y - 0.5 * (R - 1)
657  ring = max(abs(dx), abs(dy))
658  r = 2.0 * ring - 1.0
659  if r < 0:
660  if root == 5:
661  return 0
662  else:
663  return 6 * R ** 2 - 1
664  if dy == -ring:
665  i = int(r * r + ring - dx)
666  elif dx == -ring:
667  i = int(r * r + 3.0 * ring + dy)
668  elif dy == ring:
669  i = int(r * r + 5.0 * ring + dx)
670  else: # dx == ring
671  i = int(r * r + 7.0 * ring - dy)
672  if root == 5:
673  # south pole - done
674  return i
675  # north pole: use shrinking concentric squares and wind
676  # in the opposite direction
677  return 6 * R ** 2 - i - 1
678 
679  def coords(self, pixelId):
680  """Maps from a sky-pixel id to a root pixel number and x, y
681  pixel coordinates.
682  """
683  if not isinstance(pixelId, numbers.Integral):
684  raise TypeError(
685  'Sky-pixel id must be an integer')
686  R = self.resolution
687  R2 = R ** 2
688  if pixelId < 0 or pixelId >= 6 * R2:
689  raise RuntimeError(
690  'Invalid sky-pixel id %d - value must be in range [0, %d)' %
691  (pixelId, 6 * R2))
692  r = pixelId//R2
693  if r > 0 and r < 5:
694  # equatorial pixel case
695  pixelId -= R2
696  y = pixelId//(4 * R)
697  pixelId -= 4 * R * y
698  root = 1 + pixelId//R
699  x = pixelId - (root - 1) * R
700  return (root, x, y)
701  # polar root pixels
702  if r == 5:
703  root = 0
704  i = 6 * R2 - pixelId - 1
705  else:
706  root = 5
707  i = pixelId
708  s = int(math.floor(math.sqrt(i)))
709  if R & 1 == 0:
710  ring = 0.5 * (s + ((s & 1) ^ 1))
711  else:
712  ring = 0.5 * (s + (s & 1))
713  r = 2.0 * ring - 1.0
714  if r < 0.0:
715  return (root, R//2, R//2)
716  i -= r * r
717  if i <= 2.0 * ring:
718  dx = ring - i
719  dy = -ring
720  elif i <= 4.0 * ring:
721  dx = -ring
722  dy = i - 3.0 * ring
723  elif i <= 6.0 * ring:
724  dx = i - 5.0 * ring
725  dy = ring
726  else:
727  dx = ring
728  dy = 7.0 * ring - i
729  return (root, int(dx + 0.5 * (R - 1)), int(dy + 0.5 * (R - 1)))
730 
731  def _fiducialXPlane(self, root, ix):
732  assert isinstance(ix, numbers.Integral) and ix >= 0 and ix <= self.resolution
733  x = 2.0 * float(ix) / float(self.resolution) - 1.0
734  c, b = self.center[root], self.x[root]
735  v = (c[0] + x * b[0], c[1] + x * b[1], c[2] + x * b[2])
736  return geom.normalize(self.xrot[root](v, 1.0, 0.0))
737 
738  def _fiducialYPlane(self, root, iy):
739  assert isinstance(iy, numbers.Integral) and iy >= 0 and iy <= self.resolution
740  y = 2.0 * float(iy) / float(self.resolution) - 1.0
741  c, b = self.center[root], self.y[root]
742  v = (c[0] + y * b[0], c[1] + y * b[1], c[2] + y * b[2])
743  return geom.normalize(self.yrot[root](v, 1.0, 0.0))
744 
745  def _intersect(self, pixels, poly, root, box):
746  dx = box[1] - box[0]
747  dy = box[3] - box[2]
748  if dx == 1 and dy == 1:
749  i = self.id(root, box[0], box[2])
750  if poly.intersects(self.getGeometry(i)):
751  pixels.append(i)
752  return
753  if dx >= dy:
754  # Split x range and recurse
755  xsplit = box[0] + dx//2
756  p = poly.clip(self.xplane[root][xsplit][1])
757  if p is not None:
758  self._intersect(pixels, p, root, (box[0], xsplit, box[2], box[3]))
759  p = poly.clip(self.xplane[root][xsplit][2])
760  if p is not None:
761  self._intersect(pixels, p, root, (xsplit, box[1], box[2], box[3]))
762  else:
763  # Split y range and recurse
764  ysplit = box[2] + dy//2
765  p = poly.clip(self.yplane[root][ysplit][1])
766  if p is not None:
767  self._intersect(pixels, p, root, (box[0], box[1], box[2], ysplit))
768  p = poly.clip(self.yplane[root][ysplit][2])
769  if p is not None:
770  self._intersect(pixels, p, root, (box[0], box[1], ysplit, box[3]))
771 
772 
773 def createQuadSpherePixelization(policy=None):
774  policyFile = pexPolicy.DefaultPolicyFile(
775  "skypix", "QuadSpherePixelizationDictionary.paf", "policy")
776  defaults = pexPolicy.Policy.createPolicy(
777  policyFile, policyFile.getRepositoryPath())
778  if policy is None:
779  policy = pexPolicy.Policy()
780  policy.mergeDefaults(defaults)
781  # Obtain resolution and padding from policy
782  resolution = policy.get('resolutionPix')
783  padding = math.radians(policy.get('paddingArcsec') / 3600.0)
784  return QuadSpherePixelization(resolution, padding)
785 
786 
787 def imageToPolygon(wcs, widthPix, heightPix, padRad=0.0):
788  """Computes and returns a spherical convex polygon approximation to the
789  region of the unit sphere covered by an image specified with a WCS and
790  a width/height (pixels). The polygon is computed by connecting the
791  world coordinates of the 4 image corners with great circles, and can
792  optionally be padded by padRad radians.
793  """
794  # Compute image corners
795  cd = wcs.getCdMatrix()
796  xpad = math.degrees(padRad) / math.sqrt(cd[0, 0]**2 + cd[0, 1]**2)
797  ypad = math.degrees(padRad) / math.sqrt(cd[1, 0]**2 + cd[1, 1]**2)
798  xmin, ymin = -0.5 - xpad, -0.5 - ypad
799  xmax, ymax = widthPix + xpad - 0.5, heightPix + ypad - 0.5
800  # Produce a lsst.afw.coord.coordLib.Coord object for each vertex
801  coords = [wcs.pixelToSky(xmin, ymin), wcs.pixelToSky(xmax, ymin),
802  wcs.pixelToSky(xmax, ymax), wcs.pixelToSky(xmin, ymax)]
803  # Map these to python tuples containing cartesian unit vectors
804  verts = []
805  for c in coords:
806  verts.append(tuple(c.getVector()))
807  # and turn them into a spherical convex polygon
808  convex, cc = geom.convex(verts)
809  if not convex:
810  raise RuntimeError('Image corners do not form a convex polygon: ' + cc)
811  elif not cc:
812  verts.reverse()
813  return geom.SphericalConvexPolygon(verts)
lsst::afw::geom::all
bool all(CoordinateExpr< N > const &expr)
Return true if all elements are true.
Definition: CoordinateExpr.h:85
lsst::pex::policy::Policy
a container for holding hierarchical configuration data in memory.
Definition: Policy.h:169
lsst::pex::policy::DefaultPolicyFile
a representation of a default Policy file that is stored as a file in the installation directory of a...
Definition: DefaultPolicyFile.h:69
lsst.skypix.quadsphere.createQuadSpherePixelization
def createQuadSpherePixelization(policy=None)
Definition: quadsphere.py:773
lsst.skypix.quadsphere._rotY
def _rotY(v, sin, cos)
Definition: quadsphere.py:39
ast::append
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
lsst.skypix.quadsphere.QuadSpherePixelization._intersect
def _intersect(self, pixels, poly, root, box)
Definition: quadsphere.py:745
lsst.skypix.quadsphere._rotNX
def _rotNX(v, sin, cos)
Definition: quadsphere.py:47
max
int max
Definition: BoundedField.cc:99
lsst.skypix.quadsphere._rotZ
def _rotZ(v, sin, cos)
Definition: quadsphere.py:43
lsst.skypix.quadsphere.QuadSpherePixelization.__init__
def __init__(self, resolution, paddingRad)
Definition: quadsphere.py:89
lsst.skypix.quadsphere._rotNY
def _rotNY(v, sin, cos)
Definition: quadsphere.py:51
lsst.skypix.quadsphere.imageToPolygon
def imageToPolygon(wcs, widthPix, heightPix, padRad=0.0)
Definition: quadsphere.py:787
lsst.skypix.quadsphere.QuadSpherePixelization.getGeometry
def getGeometry(self, pixelId, fiducial=False)
Definition: quadsphere.py:387
lsst.skypix.quadsphere._rotX
def _rotX(v, sin, cos)
Definition: quadsphere.py:35
lsst.skypix.quadsphere._rotNZ
def _rotNZ(v, sin, cos)
Definition: quadsphere.py:55
Generated on Sat Apr 7 2018 06:23:27 for LSSTApplications by   doxygen 1.8.13