Searching...
No Matches
ringsSkyMap.py
Go to the documentation of this file.
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
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
25import struct
26import math
27import numpy as np
28
29from lsst.pex.config import Field
30import lsst.geom as geom
31from .cachingSkyMap import CachingSkyMap
32from .tractInfo import ExplicitTractInfo
33
34
35class 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
82 # calculating the ring size.
83 self._ringSize = math.pi / (config.numRings + 1) # Size of a ring in Declination (radians)
84 self._ringNums = [] # Number of tracts for each ring
85 for i in range(config.numRings):
86 startDec = self._ringSize*(i + 0.5) - 0.5*math.pi
87 stopDec = startDec + self._ringSize
88 dec = min(math.fabs(startDec), math.fabs(stopDec)) # Declination for determining division in RA
89 self._ringNums.append(int(2*math.pi*math.cos(dec)/self._ringSize) + 1)
90 numTracts = sum(self._ringNums) + 2
91 super(RingsSkyMap, self).__init__(numTracts, config, version)
92 self._raStart = self.configconfig.raStart*geom.degrees
93
94 def getRingIndices(self, index):
95 """Calculate ring indices given a numerical index of a tract.
96
97 The ring indices are the ring number and the tract number within
98 the ring.
99
100 The ring number is -1 for the south polar cap and increases to the
101 north. The north polar cap has ring number = numRings. The tract
102 number is zero for either of the polar caps.
103 """
104 if index == 0: # South polar cap
105 return -1, 0
106 if index == self._numTracts - 1: # North polar cap
107 return self.configconfig.numRings, 0
108 if index < 0 or index >= self._numTracts:
109 raise IndexError("Tract index %d is out of range [0, %d]" % (index, len(self) - 1))
110 ring = 0 # Ring number
111 tractNum = index - 1 # Tract number within ring
113 # Maintain the off-by-one bug in version=0 (DM-14809).
114 # This means that the first tract in the first ring is duplicated
115 # and the first tract in the last ring is missing.
116 while ring < self.configconfig.numRings and tractNum > self._ringNums[ring]:
117 tractNum -= self._ringNums[ring]
118 ring += 1
119 else:
120 while ring < self.configconfig.numRings and tractNum >= self._ringNums[ring]:
121 tractNum -= self._ringNums[ring]
122 ring += 1
123 return ring, tractNum
124
125 def generateTract(self, index):
126 """Generate TractInfo for the specified tract index."""
127 ringNum, tractNum = self.getRingIndices(index)
128 if ringNum == -1: # South polar cap
129 ra, dec = 0, -0.5*math.pi
130 elif ringNum == self.configconfig.numRings: # North polar cap
131 ra, dec = 0, 0.5*math.pi
132 else:
133 dec = self._ringSize*(ringNum + 1) - 0.5*math.pi
136
137 center = geom.SpherePoint(ra, dec, geom.radians)
138 wcs = self._wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=center)
141 wcs)
142
143 def _decToRingNum(self, dec):
144 """Calculate ring number from Declination.
145
146 Parameters
147 ----------
148 dec : lsst.geom.Angle
149 Declination.
150
151 Returns
152 -------
153 ringNum : int
154 Ring number: -1 for the south polar cap, and increasing to the
155 north, ending with numRings for the north polar cap.
156 """
157 firstRingStart = self._ringSize*0.5 - 0.5*math.pi
158 if dec < firstRingStart:
159 # Southern cap
160 return -1
161 elif dec > firstRingStart*-1:
162 # Northern cap
163 return self.configconfig.numRings
165
166 def _raToTractNum(self, ra, ringNum):
167 """Calculate tract number from the Right Ascension.
168
169 Parameters
170 ----------
171 ra : lsst.geom.Angle
172 Right Ascension.
173 ringNum : int
174 Ring number (from _decToRingNum).
175
176 Returns
177 -------
178 tractNum : int
179 Tract number within the ring (starts at 0 for the tract at
180 raStart).
181 """
182 if ringNum in (-1, self.configconfig.numRings):
183 return 0
184 assert ringNum in range(self.configconfig.numRings)
185 tractNum = int((ra - self._raStart).wrap().asRadians()
186 / (2*math.pi/self._ringNums[ringNum]) + 0.5)
187 return 0 if tractNum == self._ringNums[ringNum] else tractNum # Allow wraparound
188
189 def findTract(self, coord):
190 ringNum = self._decToRingNum(coord.getLatitude())
191 if ringNum == -1:
192 # Southern cap
193 return self[0]
194 if ringNum == self.configconfig.numRings:
195 # Northern cap
196 return self[self._numTracts - 1]
197 tractNum = self._raToTractNum(coord.getLongitude(), ringNum)
198
199 if self._version_version_version == 0 and tractNum == 0 and ringNum != 0:
200 # Account for off-by-one error in getRingIndices
201 # Note that this means that tract 1 gets duplicated.
202 ringNum += 1
203
204 index = sum(self._ringNums[:ringNum], tractNum + 1) # Allow 1 for south pole
205 return self[index]
206
207 def findTractIdArray(self, ra, dec, degrees=False):
208 if degrees:
211 else:
212 _ra = np.atleast_1d(ra)
213 _dec = np.atleast_1d(dec)
214
215 # Set default to -1 to distinguish from polar tracts
216 indexes = np.full(_ra.size, -1, dtype=np.int32)
217
218 # Do the dec search
219 firstRingStart = self._ringSize*0.5 - 0.5*np.pi
220 ringNums = np.zeros(len(_dec), dtype=np.int32)
221
222 ringNums[_dec < firstRingStart] = -1
223 ringNums[_dec > -1*firstRingStart] = self.configconfig.numRings
224
225 mid = (_dec >= firstRingStart) & (_dec <= -1*firstRingStart)
226 ringNums[mid] = ((_dec[mid] - firstRingStart)/self._ringSize).astype(np.int32)
227
228 indexes[ringNums == -1] = 0
229 indexes[ringNums == self.configconfig.numRings] = self._numTracts - 1
230
231 # We now do the full lookup for all non-polar tracts that have not
232 # been set.
233 inRange, = np.where(indexes < 0)
234
235 # Do the ra search
236 _ringNumArray = np.array(self._ringNums)
237 _ringCumulative = np.cumsum(np.insert(_ringNumArray, 0, 0))
238
239 deltaWrap = (_ra[inRange] - self._raStart.asRadians()) % (2.*np.pi)
240 tractNum = (deltaWrap/(2.*np.pi/_ringNumArray[ringNums[inRange]]) + 0.5).astype(np.int32)
241 # Allow wraparound
242 tractNum[tractNum == _ringNumArray[ringNums[inRange]]] = 0
243
244 if self._version_version_version == 0:
245 # Account for off-by-one error in getRingIndices
246 # Note that this means tract 1 gets duplicated
247 offByOne, = np.where((tractNum == 0)
248 & (ringNums[inRange] != 0))
249 ringNums[inRange[offByOne]] += 1
250
251 indexes[inRange] = _ringCumulative[ringNums[inRange]] + tractNum + 1
252
253 return indexes
254
255 def findAllTracts(self, coord):
256 """Find all tracts which include the specified coord.
257
258 Parameters
259 ----------
260 coord : lsst.geom.SpherePoint
261 ICRS sky coordinate to search for.
262
263 Returns
264 -------
265 tractList : list of TractInfo
266 The tracts which include the specified coord.
267 """
268 ringNum = self._decToRingNum(coord.getLatitude())
269
270 tractList = list()
271 # ringNum denotes the closest ring to the specified coord
272 # I will check adjacent rings which may include the specified coord
273 for r in [ringNum - 1, ringNum, ringNum + 1]:
274 if r < 0 or r >= self.configconfig.numRings:
275 # Poles will be checked explicitly outside this loop
276 continue
277 tractNum = self._raToTractNum(coord.getLongitude(), r)
278 # Adjacent tracts will also be checked.
279 for t in [tractNum - 1, tractNum, tractNum + 1]:
280 # Wrap over raStart
281 if t < 0:
282 t = t + self._ringNums[r]
283 elif t > self._ringNums[r] - 1:
284 t = t - self._ringNums[r]
285
286 extra = 0
287 if self._version_version_version == 0 and t == 0 and r != 0:
288 # Account for off-by-one error in getRingIndices
289 # Note that this means that tract 1 gets duplicated.
290 extra = 1
291
292 index = sum(self._ringNums[:r + extra], t + 1) # Allow 1 for south pole
293 tract = self[index]
294 if tract.contains(coord):
295 tractList.append(tract)
296
297 # Always check tracts at poles
298 # Southern cap is 0, Northern cap is the last entry in self
299 for entry in [0, len(self)-1]:
300 tract = self[entry]
301 if tract.contains(coord):
302 tractList.append(tract)
303
304 return tractList
305
306 def findTractPatchList(self, coordList):
307 retList = []
308 for coord in coordList:
309 for tractInfo in self.findAllTracts(coord):
310 patchList = tractInfo.findPatchList(coordList)
311 if patchList and not (tractInfo, patchList) in retList:
312 retList.append((tractInfo, patchList))
313 return retList
314