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
ringsSkyMap.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010, 2012 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__ = ["RingsSkyMapConfig", "RingsSkyMap"]
24 
25 import struct
26 import math
27 import numpy as np
28 
29 from lsst.pex.config import Field
30 import lsst.geom as geom
31 from .cachingSkyMap import CachingSkyMap
32 from .tractInfo import ExplicitTractInfo
33 
34 
35 class RingsSkyMapConfig(CachingSkyMap.ConfigClass):
36  """Configuration for the RingsSkyMap"""
37  numRings = Field(dtype=int, doc="Number of rings", check=lambda x: x > 0)
38  raStart = Field(dtype=float, default=0.0, doc="Starting center RA for each ring (degrees)",
39  check=lambda x: x >= 0.0 and x < 360.0)
40 
41 
43  """Rings sky map pixelization.
44 
45  We divide the sphere into N rings of Declination, plus the two polar
46  caps, which sets the size of the individual tracts. The rings are
47  divided in RA into an integral number of tracts of this size; this
48  division is made at the Declination closest to zero so as to ensure
49  full overlap.
50 
51  Rings are numbered in the rings from south to north. The south pole cap is
52  ``tract=0``, then the tract at ``raStart`` in the southernmost ring is
53  ``tract=1``. Numbering continues (in the positive RA direction) around that
54  ring and then continues in the same fashion with the next ring north, and
55  so on until all reaching the north pole cap, which is
56  ``tract=len(skymap) - 1``.
57 
58  However, ``version=0`` had a bug in the numbering of the tracts: the first
59  and last tracts in the first (southernmost) ring were identical, and the
60  first tract in the last (northernmost) ring was missing. When using
61  ``version=0``, these tracts remain missing in order to preserve the
62  numbering scheme.
63 
64  Parameters
65  ----------
66  config : `lsst.skymap.RingsSkyMapConfig`
67  The configuration for this SkyMap.
68  version : `int`, optional
69  Software version of this class, to retain compatibility with old
70  verisons. ``version=0`` covers the period from first implementation
71  until DM-14809, at which point bugs were identified in the numbering
72  of tracts (affecting only tracts at RA=0). ``version=1`` uses the
73  post-DM-14809 tract numbering.
74  """
75  ConfigClass = RingsSkyMapConfig
76  _version = (1, 0) # for pickle
77 
78  def __init__(self, config, version=1):
79  assert version in (0, 1), "Unrecognised version: %s" % (version,)
80  # We count rings from south to north
81  # Note: pole caps together count for one additional ring when calculating the ring size
82  self._ringSize_ringSize = math.pi / (config.numRings + 1) # Size of a ring in Declination (radians)
83  self._ringNums_ringNums = [] # Number of tracts for each ring
84  for i in range(config.numRings):
85  startDec = self._ringSize_ringSize*(i + 0.5) - 0.5*math.pi
86  stopDec = startDec + self._ringSize_ringSize
87  dec = min(math.fabs(startDec), math.fabs(stopDec)) # Declination for determining division in RA
88  self._ringNums_ringNums.append(int(2*math.pi*math.cos(dec)/self._ringSize_ringSize) + 1)
89  numTracts = sum(self._ringNums_ringNums) + 2
90  super(RingsSkyMap, self).__init__(numTracts, config, version)
91  self._raStart_raStart = self.configconfig.raStart*geom.degrees
92 
93  def getRingIndices(self, index):
94  """Calculate ring indices given a numerical index of a tract.
95 
96  The ring indices are the ring number and the tract number within
97  the ring.
98 
99  The ring number is -1 for the south polar cap and increases to the
100  north. The north polar cap has ring number = numRings. The tract
101  number is zero for either of the polar caps.
102  """
103  if index == 0: # South polar cap
104  return -1, 0
105  if index == self._numTracts_numTracts - 1: # North polar cap
106  return self.configconfig.numRings, 0
107  if index < 0 or index >= self._numTracts_numTracts:
108  raise IndexError("Tract index %d is out of range [0, %d]" % (index, len(self) - 1))
109  ring = 0 # Ring number
110  tractNum = index - 1 # Tract number within ring
111  if self._version_version_version_version == 0:
112  # Maintain the off-by-one bug in version=0 (DM-14809).
113  # This means that the first tract in the first ring is duplicated
114  # and the first tract in the last ring is missing.
115  while ring < self.configconfig.numRings and tractNum > self._ringNums_ringNums[ring]:
116  tractNum -= self._ringNums_ringNums[ring]
117  ring += 1
118  else:
119  while ring < self.configconfig.numRings and tractNum >= self._ringNums_ringNums[ring]:
120  tractNum -= self._ringNums_ringNums[ring]
121  ring += 1
122  return ring, tractNum
123 
124  def generateTract(self, index):
125  """Generate TractInfo for the specified tract index."""
126  ringNum, tractNum = self.getRingIndicesgetRingIndices(index)
127  if ringNum == -1: # South polar cap
128  ra, dec = 0, -0.5*math.pi
129  elif ringNum == self.configconfig.numRings: # North polar cap
130  ra, dec = 0, 0.5*math.pi
131  else:
132  dec = self._ringSize_ringSize*(ringNum + 1) - 0.5*math.pi
133  ra = ((2*math.pi*tractNum/self._ringNums_ringNums[ringNum])*geom.radians
134  + self._raStart_raStart).wrap().asRadians()
135 
136  center = geom.SpherePoint(ra, dec, geom.radians)
137  wcs = self._wcsFactory_wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=center)
138  return ExplicitTractInfo(index, self._tractBuilder_tractBuilder, center,
139  0.5*self._ringSize_ringSize*geom.radians, self.configconfig.tractOverlap*geom.degrees,
140  wcs)
141 
142  def _decToRingNum(self, dec):
143  """Calculate ring number from Declination.
144 
145  Parameters
146  ----------
147  dec : `lsst.geom.Angle`
148  Declination.
149 
150  Returns
151  -------
152  ringNum : `int`
153  Ring number: -1 for the south polar cap, and increasing to the
154  north, ending with ``numRings`` for the north polar cap.
155  """
156  firstRingStart = self._ringSize_ringSize*0.5 - 0.5*math.pi
157  if dec < firstRingStart:
158  # Southern cap
159  return -1
160  elif dec > firstRingStart*-1:
161  # Northern cap
162  return self.configconfig.numRings
163  return int((dec.asRadians() - firstRingStart)/self._ringSize_ringSize)
164 
165  def _raToTractNum(self, ra, ringNum):
166  """Calculate tract number from the Right Ascension.
167 
168  Parameters
169  ----------
170  ra : `lsst.geom.Angle`
171  Right Ascension.
172  ringNum : `int`
173  Ring number (from ``_decToRingNum``).
174 
175  Returns
176  -------
177  tractNum : `int`
178  Tract number within the ring (starts at 0 for the tract at raStart).
179  """
180  if ringNum in (-1, self.configconfig.numRings):
181  return 0
182  assert ringNum in range(self.configconfig.numRings)
183  tractNum = int((ra - self._raStart_raStart).wrap().asRadians()
184  / (2*math.pi/self._ringNums_ringNums[ringNum]) + 0.5)
185  return 0 if tractNum == self._ringNums_ringNums[ringNum] else tractNum # Allow wraparound
186 
187  def findTract(self, coord):
188  ringNum = self._decToRingNum_decToRingNum(coord.getLatitude())
189  if ringNum == -1:
190  # Southern cap
191  return self[0]
192  if ringNum == self.configconfig.numRings:
193  # Northern cap
194  return self[self._numTracts_numTracts - 1]
195  tractNum = self._raToTractNum_raToTractNum(coord.getLongitude(), ringNum)
196 
197  if self._version_version_version_version == 0 and tractNum == 0 and ringNum != 0:
198  # Account for off-by-one error in getRingIndices
199  # Note that this means that tract 1 gets duplicated.
200  ringNum += 1
201 
202  index = sum(self._ringNums_ringNums[:ringNum], tractNum + 1) # Allow 1 for south pole
203  return self[index]
204 
205  def findTractIdArray(self, ra, dec, degrees=False):
206  if degrees:
207  _ra = np.deg2rad(ra)
208  _dec = np.deg2rad(dec)
209  else:
210  _ra = np.atleast_1d(ra)
211  _dec = np.atleast_1d(dec)
212 
213  # Set default to -1 to distinguish from polar tracts
214  indexes = np.full(_ra.size, -1, dtype=np.int32)
215 
216  # Do the dec search
217  firstRingStart = self._ringSize_ringSize*0.5 - 0.5*np.pi
218  ringNums = np.zeros(len(_dec), dtype=np.int32)
219 
220  ringNums[_dec < firstRingStart] = -1
221  ringNums[_dec > -1*firstRingStart] = self.configconfig.numRings
222 
223  mid = (_dec >= firstRingStart) & (_dec <= -1*firstRingStart)
224  ringNums[mid] = ((_dec[mid] - firstRingStart)/self._ringSize_ringSize).astype(np.int32)
225 
226  indexes[ringNums == -1] = 0
227  indexes[ringNums == self.configconfig.numRings] = self._numTracts_numTracts - 1
228 
229  # We now do the full lookup for all non-polar tracts that have not been set.
230  inRange, = np.where(indexes < 0)
231 
232  # Do the ra search
233  _ringNumArray = np.array(self._ringNums_ringNums)
234  _ringCumulative = np.cumsum(np.insert(_ringNumArray, 0, 0))
235 
236  deltaWrap = (_ra[inRange] - self._raStart_raStart.asRadians()) % (2.*np.pi)
237  tractNum = (deltaWrap/(2.*np.pi/_ringNumArray[ringNums[inRange]]) + 0.5).astype(np.int32)
238  # Allow wraparound
239  tractNum[tractNum == _ringNumArray[ringNums[inRange]]] = 0
240 
241  if self._version_version_version_version == 0:
242  # Account for off-by-one error in getRingIndices
243  # Note that this means tract 1 gets duplicated
244  offByOne, = np.where((tractNum == 0)
245  & (ringNums[inRange] != 0))
246  ringNums[inRange[offByOne]] += 1
247 
248  indexes[inRange] = _ringCumulative[ringNums[inRange]] + tractNum + 1
249 
250  return indexes
251 
252  def findAllTracts(self, coord):
253  """Find all tracts which include the specified coord.
254 
255  Parameters
256  ----------
257  coord : `lsst.geom.SpherePoint`
258  ICRS sky coordinate to search for.
259 
260  Returns
261  -------
262  tractList : `list` of `TractInfo`
263  The tracts which include the specified coord.
264  """
265  ringNum = self._decToRingNum_decToRingNum(coord.getLatitude())
266 
267  tractList = list()
268  # ringNum denotes the closest ring to the specified coord
269  # I will check adjacent rings which may include the specified coord
270  for r in [ringNum - 1, ringNum, ringNum + 1]:
271  if r < 0 or r >= self.configconfig.numRings:
272  # Poles will be checked explicitly outside this loop
273  continue
274  tractNum = self._raToTractNum_raToTractNum(coord.getLongitude(), r)
275  # Adjacent tracts will also be checked.
276  for t in [tractNum - 1, tractNum, tractNum + 1]:
277  # Wrap over raStart
278  if t < 0:
279  t = t + self._ringNums_ringNums[r]
280  elif t > self._ringNums_ringNums[r] - 1:
281  t = t - self._ringNums_ringNums[r]
282 
283  extra = 0
284  if self._version_version_version_version == 0 and t == 0 and r != 0:
285  # Account for off-by-one error in getRingIndices
286  # Note that this means that tract 1 gets duplicated.
287  extra = 1
288 
289  index = sum(self._ringNums_ringNums[:r + extra], t + 1) # Allow 1 for south pole
290  tract = self[index]
291  if tract.contains(coord):
292  tractList.append(tract)
293 
294  # Always check tracts at poles
295  # Southern cap is 0, Northern cap is the last entry in self
296  for entry in [0, len(self)-1]:
297  tract = self[entry]
298  if tract.contains(coord):
299  tractList.append(tract)
300 
301  return tractList
302 
303  def findTractPatchList(self, coordList):
304  retList = []
305  for coord in coordList:
306  for tractInfo in self.findAllTractsfindAllTracts(coord):
307  patchList = tractInfo.findPatchList(coordList)
308  if patchList and not (tractInfo, patchList) in retList:
309  retList.append((tractInfo, patchList))
310  return retList
311 
312  def updateSha1(self, sha1):
313  """Add subclass-specific state or configuration options to the SHA1."""
314  sha1.update(struct.pack("<id", self.configconfig.numRings, self.configconfig.raStart))
int min
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
def _raToTractNum(self, ra, ringNum)
Definition: ringsSkyMap.py:165
def __init__(self, config, version=1)
Definition: ringsSkyMap.py:78
def findTractIdArray(self, ra, dec, degrees=False)
Definition: ringsSkyMap.py:205
def findTractPatchList(self, coordList)
Definition: ringsSkyMap.py:303
daf::base::PropertyList * list
Definition: fits.cc:913
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
std::shared_ptr< afw::geom::SkyWcs > makeWcs(SipForwardTransform const &sipForward, SipReverseTransform const &sipReverse, geom::SpherePoint const &skyOrigin)
Create a new TAN SIP Wcs from a pair of SIP transforms and the sky origin.
def wrap(ctrl)
Definition: wrap.py:302