23 __all__ = [
"RingsSkyMapConfig",
"RingsSkyMap"]
31 from .cachingSkyMap
import CachingSkyMap
32 from .tractInfo
import ExplicitTractInfo
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)
43 """Rings sky map pixelization.
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
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``.
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
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.
75 ConfigClass = RingsSkyMapConfig
79 assert version
in (0, 1),
"Unrecognised version: %s" % (version,)
82 self.
_ringSize_ringSize = math.pi / (config.numRings + 1)
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))
89 numTracts = sum(self.
_ringNums_ringNums) + 2
90 super(RingsSkyMap, self).
__init__(numTracts, config, version)
94 """Calculate ring indices given a numerical index of a tract.
96 The ring indices are the ring number and the tract number within
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.
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))
115 while ring < self.
configconfig.numRings
and tractNum > self.
_ringNums_ringNums[ring]:
116 tractNum -= self.
_ringNums_ringNums[ring]
119 while ring < self.
configconfig.numRings
and tractNum >= self.
_ringNums_ringNums[ring]:
120 tractNum -= self.
_ringNums_ringNums[ring]
122 return ring, tractNum
125 """Generate TractInfo for the specified tract index."""
128 ra, dec = 0, -0.5*math.pi
129 elif ringNum == self.
configconfig.numRings:
130 ra, dec = 0, 0.5*math.pi
132 dec = self.
_ringSize_ringSize*(ringNum + 1) - 0.5*math.pi
133 ra = ((2*math.pi*tractNum/self.
_ringNums_ringNums[ringNum])*geom.radians
139 0.5*self.
_ringSize_ringSize*geom.radians, self.
configconfig.tractOverlap*geom.degrees,
142 def _decToRingNum(self, dec):
143 """Calculate ring number from Declination.
147 dec : `lsst.geom.Angle`
153 Ring number: -1 for the south polar cap, and increasing to the
154 north, ending with ``numRings`` for the north polar cap.
156 firstRingStart = self.
_ringSize_ringSize*0.5 - 0.5*math.pi
157 if dec < firstRingStart:
160 elif dec > firstRingStart*-1:
162 return self.
configconfig.numRings
163 return int((dec.asRadians() - firstRingStart)/self.
_ringSize_ringSize)
165 def _raToTractNum(self, ra, ringNum):
166 """Calculate tract number from the Right Ascension.
170 ra : `lsst.geom.Angle`
173 Ring number (from ``_decToRingNum``).
178 Tract number within the ring (starts at 0 for the tract at raStart).
180 if ringNum
in (-1, self.
configconfig.numRings):
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
188 ringNum = self.
_decToRingNum_decToRingNum(coord.getLatitude())
192 if ringNum == self.
configconfig.numRings:
195 tractNum = self.
_raToTractNum_raToTractNum(coord.getLongitude(), ringNum)
202 index = sum(self.
_ringNums_ringNums[:ringNum], tractNum + 1)
208 _dec = np.deg2rad(dec)
210 _ra = np.atleast_1d(ra)
211 _dec = np.atleast_1d(dec)
214 indexes = np.full(_ra.size, -1, dtype=np.int32)
217 firstRingStart = self.
_ringSize_ringSize*0.5 - 0.5*np.pi
218 ringNums = np.zeros(len(_dec), dtype=np.int32)
220 ringNums[_dec < firstRingStart] = -1
221 ringNums[_dec > -1*firstRingStart] = self.
configconfig.numRings
223 mid = (_dec >= firstRingStart) & (_dec <= -1*firstRingStart)
224 ringNums[mid] = ((_dec[mid] - firstRingStart)/self.
_ringSize_ringSize).astype(np.int32)
226 indexes[ringNums == -1] = 0
227 indexes[ringNums == self.
configconfig.numRings] = self.
_numTracts_numTracts - 1
230 inRange, = np.where(indexes < 0)
233 _ringNumArray = np.array(self.
_ringNums_ringNums)
234 _ringCumulative = np.cumsum(np.insert(_ringNumArray, 0, 0))
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)
239 tractNum[tractNum == _ringNumArray[ringNums[inRange]]] = 0
244 offByOne, = np.where((tractNum == 0)
245 & (ringNums[inRange] != 0))
246 ringNums[inRange[offByOne]] += 1
248 indexes[inRange] = _ringCumulative[ringNums[inRange]] + tractNum + 1
253 """Find all tracts which include the specified coord.
257 coord : `lsst.geom.SpherePoint`
258 ICRS sky coordinate to search for.
262 tractList : `list` of `TractInfo`
263 The tracts which include the specified coord.
265 ringNum = self.
_decToRingNum_decToRingNum(coord.getLatitude())
270 for r
in [ringNum - 1, ringNum, ringNum + 1]:
271 if r < 0
or r >= self.
configconfig.numRings:
274 tractNum = self.
_raToTractNum_raToTractNum(coord.getLongitude(), r)
276 for t
in [tractNum - 1, tractNum, tractNum + 1]:
289 index = sum(self.
_ringNums_ringNums[:r + extra], t + 1)
291 if tract.contains(coord):
292 tractList.append(tract)
296 for entry
in [0, len(self)-1]:
298 if tract.contains(coord):
299 tractList.append(tract)
305 for coord
in coordList:
307 patchList = tractInfo.findPatchList(coordList)
308 if patchList
and not (tractInfo, patchList)
in retList:
309 retList.append((tractInfo, patchList))
313 """Add subclass-specific state or configuration options to the SHA1."""
314 sha1.update(struct.pack(
"<id", self.
configconfig.numRings, self.
configconfig.raStart))
Point in an unspecified spherical coordinate system.
def getRingIndices(self, index)
def _raToTractNum(self, ra, ringNum)
def _decToRingNum(self, dec)
def __init__(self, config, version=1)
def findTractIdArray(self, ra, dec, degrees=False)
def updateSha1(self, sha1)
def findTract(self, coord)
def findTractPatchList(self, coordList)
def generateTract(self, index)
def findAllTracts(self, coord)
daf::base::PropertyList * list
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
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.