Coverage for python / lsst / skymap / ringsSkyMap.py: 12%

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 

25import struct 

26import math 

27import numpy as np 

28 

29from lsst.pex.config import Field 

30import lsst.geom as geom 

31import lsst.sphgeom as sphgeom 

32from .cachingSkyMap import CachingSkyMap 

33from .tractInfo import ExplicitTractInfo 

34 

35 

36class RingsSkyMapConfig(CachingSkyMap.ConfigClass): 

37 """Configuration for the RingsSkyMap""" 

38 numRings = Field(dtype=int, doc="Number of rings", check=lambda x: x > 0) 

39 raStart = Field(dtype=float, default=0.0, doc="Starting center RA for each ring (degrees)", 

40 check=lambda x: x >= 0.0 and x < 360.0) 

41 

42 

43class RingsSkyMap(CachingSkyMap): 

44 """Rings sky map pixelization. 

45 

46 We divide the sphere into N rings of Declination, plus the two polar 

47 caps, which sets the size of the individual tracts. The rings are 

48 divided in RA into an integral number of tracts of this size; this 

49 division is made at the Declination closest to zero so as to ensure 

50 full overlap. 

51 

52 Rings are numbered in the rings from south to north. The south pole cap is 

53 ``tract=0``, then the tract at ``raStart`` in the southernmost ring is 

54 ``tract=1``. Numbering continues (in the positive RA direction) around that 

55 ring and then continues in the same fashion with the next ring north, and 

56 so on until all reaching the north pole cap, which is 

57 ``tract=len(skymap) - 1``. 

58 

59 However, ``version=0`` had a bug in the numbering of the tracts: the first 

60 and last tracts in the first (southernmost) ring were identical, and the 

61 first tract in the last (northernmost) ring was missing. When using 

62 ``version=0``, these tracts remain missing in order to preserve the 

63 numbering scheme. 

64 

65 Parameters 

66 ---------- 

67 config : `lsst.skymap.RingsSkyMapConfig` 

68 The configuration for this SkyMap. 

69 version : `int`, optional 

70 Software version of this class, to retain compatibility with old 

71 verisons. ``version=0`` covers the period from first implementation 

72 until DM-14809, at which point bugs were identified in the numbering 

73 of tracts (affecting only tracts at RA=0). ``version=1`` uses the 

74 post-DM-14809 tract numbering. 

75 """ 

76 ConfigClass = RingsSkyMapConfig 

77 _version = (1, 0) # for pickle 

78 

79 def __init__(self, config, version=1): 

80 assert version in (0, 1), "Unrecognised version: %s" % (version,) 

81 # We count rings from south to north 

82 # Note: pole caps together count for one additional ring when 

83 # calculating the ring size. 

84 self._ringSize = math.pi / (config.numRings + 1) # Size of a ring in Declination (radians) 

85 self._ringNums = [] # Number of tracts for each ring 

86 for i in range(config.numRings): 

87 startDec = self._ringSize*(i + 0.5) - 0.5*math.pi 

88 stopDec = startDec + self._ringSize 

89 dec = min(math.fabs(startDec), math.fabs(stopDec)) # Declination for determining division in RA 

90 self._ringNums.append(int(2*math.pi*math.cos(dec)/self._ringSize) + 1) 

91 numTracts = sum(self._ringNums) + 2 

92 super(RingsSkyMap, self).__init__(numTracts, config, version) 

93 self._raStart = self.config.raStart*geom.degrees 

94 

95 def getRingIndices(self, index): 

96 """Calculate ring indices given a numerical index of a tract. 

97 

98 The ring indices are the ring number and the tract number within 

99 the ring. 

100 

101 The ring number is -1 for the south polar cap and increases to the 

102 north. The north polar cap has ring number = numRings. The tract 

103 number is zero for either of the polar caps. 

104 """ 

105 if index == 0: # South polar cap 

106 return -1, 0 

107 if index == self._numTracts - 1: # North polar cap 

108 return self.config.numRings, 0 

109 if index < 0 or index >= self._numTracts: 

110 raise IndexError("Tract index %d is out of range [0, %d]" % (index, len(self) - 1)) 

111 ring = 0 # Ring number 

112 tractNum = index - 1 # Tract number within ring 

113 if self._version == 0: 

114 # Maintain the off-by-one bug in version=0 (DM-14809). 

115 # This means that the first tract in the first ring is duplicated 

116 # and the first tract in the last ring is missing. 

117 while ring < self.config.numRings and tractNum > self._ringNums[ring]: 

118 tractNum -= self._ringNums[ring] 

119 ring += 1 

120 else: 

121 while ring < self.config.numRings and tractNum >= self._ringNums[ring]: 

122 tractNum -= self._ringNums[ring] 

123 ring += 1 

124 return ring, tractNum 

125 

126 def generateTract(self, index): 

127 """Generate TractInfo for the specified tract index.""" 

128 ringNum, tractNum = self.getRingIndices(index) 

129 if ringNum == -1: # South polar cap 

130 ra, dec = 0, -0.5*math.pi 

131 elif ringNum == self.config.numRings: # North polar cap 

132 ra, dec = 0, 0.5*math.pi 

133 else: 

134 dec = self._ringSize*(ringNum + 1) - 0.5*math.pi 

135 ra = ((2*math.pi*tractNum/self._ringNums[ringNum])*geom.radians 

136 + self._raStart).wrap().asRadians() 

137 

138 center = geom.SpherePoint(ra, dec, geom.radians) 

139 wcs = self._wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=center) 

140 

141 raMin, raMax, decMin, decMax = self.getRaDecRange(index) 

142 

143 innerBoxCorners = [ 

144 sphgeom.LonLat(lon=sphgeom.NormalizedAngle(raMin), lat=sphgeom.Angle(decMin)), 

145 sphgeom.LonLat(lon=sphgeom.NormalizedAngle(raMax), lat=sphgeom.Angle(decMax)), 

146 ] 

147 

148 return ExplicitTractInfo( 

149 index, 

150 self._tractBuilder, 

151 center, 

152 0.5*self._ringSize*geom.radians, 

153 self.config.tractOverlap*geom.degrees, 

154 wcs, 

155 innerBoxCorners=innerBoxCorners, 

156 ) 

157 

158 def getRaDecRange(self, index): 

159 """Get the ra and dec ranges for the inner region of a 

160 specified tract index. 

161 

162 Parameters 

163 ---------- 

164 index : `int` 

165 Tract index number. 

166 

167 Returns 

168 ------- 

169 raMin, raMax, decMin, decMax : `lsst.geom.Angle` 

170 RA/Dec boundaries of the inner region. 

171 """ 

172 ringNum, tractNum = self.getRingIndices(index) 

173 

174 if ringNum == -1: 

175 # South polar cap. 

176 decMin = (-np.pi/2.)*geom.radians 

177 decMax = (-np.pi/2. + (ringNum + 1.5)*self._ringSize)*geom.radians 

178 raMin = 0.0*geom.radians 

179 raMax = 2.*np.pi*geom.radians 

180 elif ringNum == self.config.numRings: 

181 # North polar cap. 

182 decMin = (-np.pi/2. + (ringNum + 0.5)*self._ringSize)*geom.radians 

183 decMax = (np.pi/2.)*geom.radians 

184 raMin = 0.0*geom.radians 

185 raMax = 2.*np.pi*geom.radians 

186 else: 

187 decMin = (-np.pi/2. + (ringNum + 0.5)*self._ringSize)*geom.radians 

188 decMax = (-np.pi/2. + (ringNum + 1.5)*self._ringSize)*geom.radians 

189 

190 deltaStart = ((tractNum - 0.5)*2.*np.pi/self._ringNums[ringNum])*geom.radians 

191 deltaStop = ((tractNum + 0.5)*2.*np.pi/self._ringNums[ringNum])*geom.radians 

192 

193 raMin = (self._raStart + deltaStart).wrap() 

194 raMax = (self._raStart + deltaStop).wrap() 

195 

196 return raMin, raMax, decMin, decMax 

197 

198 def _decToRingNum(self, dec): 

199 """Calculate ring number from Declination. 

200 

201 Parameters 

202 ---------- 

203 dec : `lsst.geom.Angle` 

204 Declination. 

205 

206 Returns 

207 ------- 

208 ringNum : `int` 

209 Ring number: -1 for the south polar cap, and increasing to the 

210 north, ending with ``numRings`` for the north polar cap. 

211 """ 

212 firstRingStart = self._ringSize*0.5 - 0.5*math.pi 

213 if dec < firstRingStart: 

214 # Southern cap 

215 return -1 

216 elif dec > firstRingStart*-1: 

217 # Northern cap 

218 return self.config.numRings 

219 return int((dec.asRadians() - firstRingStart)/self._ringSize) 

220 

221 def _raToTractNum(self, ra, ringNum): 

222 """Calculate tract number from the Right Ascension. 

223 

224 Parameters 

225 ---------- 

226 ra : `lsst.geom.Angle` 

227 Right Ascension. 

228 ringNum : `int` 

229 Ring number (from ``_decToRingNum``). 

230 

231 Returns 

232 ------- 

233 tractNum : `int` 

234 Tract number within the ring (starts at 0 for the tract at 

235 ``raStart``). 

236 """ 

237 if ringNum in (-1, self.config.numRings): 

238 return 0 

239 assert ringNum in range(self.config.numRings) 

240 tractNum = int((ra - self._raStart).wrap().asRadians() 

241 / (2*math.pi/self._ringNums[ringNum]) + 0.5) 

242 return 0 if tractNum == self._ringNums[ringNum] else tractNum # Allow wraparound 

243 

244 def findTract(self, coord): 

245 ringNum = self._decToRingNum(coord.getLatitude()) 

246 if ringNum == -1: 

247 # Southern cap 

248 return self[0] 

249 if ringNum == self.config.numRings: 

250 # Northern cap 

251 return self[self._numTracts - 1] 

252 tractNum = self._raToTractNum(coord.getLongitude(), ringNum) 

253 

254 if self._version == 0 and tractNum == 0 and ringNum != 0: 

255 # Account for off-by-one error in getRingIndices 

256 # Note that this means that tract 1 gets duplicated. 

257 ringNum += 1 

258 

259 index = sum(self._ringNums[:ringNum], tractNum + 1) # Allow 1 for south pole 

260 return self[index] 

261 

262 def findTractIdArray(self, ra, dec, degrees=False): 

263 if degrees: 

264 _ra = np.deg2rad(ra) 

265 _dec = np.deg2rad(dec) 

266 else: 

267 _ra = np.atleast_1d(ra) 

268 _dec = np.atleast_1d(dec) 

269 

270 # Set default to -1 to distinguish from polar tracts 

271 indexes = np.full(_ra.size, -1, dtype=np.int32) 

272 

273 # Do the dec search 

274 firstRingStart = self._ringSize*0.5 - 0.5*np.pi 

275 ringNums = np.zeros(len(_dec), dtype=np.int32) 

276 

277 ringNums[_dec < firstRingStart] = -1 

278 ringNums[_dec > -1*firstRingStart] = self.config.numRings 

279 

280 mid = (_dec >= firstRingStart) & (_dec <= -1*firstRingStart) 

281 ringNums[mid] = ((_dec[mid] - firstRingStart)/self._ringSize).astype(np.int32) 

282 

283 indexes[ringNums == -1] = 0 

284 indexes[ringNums == self.config.numRings] = self._numTracts - 1 

285 

286 # We now do the full lookup for all non-polar tracts that have not 

287 # been set. 

288 inRange, = np.where(indexes < 0) 

289 

290 # Do the ra search 

291 _ringNumArray = np.array(self._ringNums) 

292 _ringCumulative = np.cumsum(np.insert(_ringNumArray, 0, 0)) 

293 

294 deltaWrap = (_ra[inRange] - self._raStart.asRadians()) % (2.*np.pi) 

295 tractNum = (deltaWrap/(2.*np.pi/_ringNumArray[ringNums[inRange]]) + 0.5).astype(np.int32) 

296 # Allow wraparound 

297 tractNum[tractNum == _ringNumArray[ringNums[inRange]]] = 0 

298 

299 if self._version == 0: 

300 # Account for off-by-one error in getRingIndices 

301 # Note that this means tract 1 gets duplicated 

302 offByOne, = np.where((tractNum == 0) 

303 & (ringNums[inRange] != 0)) 

304 ringNums[inRange[offByOne]] += 1 

305 

306 indexes[inRange] = _ringCumulative[ringNums[inRange]] + tractNum + 1 

307 

308 return indexes 

309 

310 def findAllTracts(self, coord): 

311 """Find all tracts which include the specified coord. 

312 

313 Parameters 

314 ---------- 

315 coord : `lsst.geom.SpherePoint` 

316 ICRS sky coordinate to search for. 

317 

318 Returns 

319 ------- 

320 tractList : `list` of `TractInfo` 

321 The tracts which include the specified coord. 

322 """ 

323 ringNum = self._decToRingNum(coord.getLatitude()) 

324 

325 tractList = list() 

326 # ringNum denotes the closest ring to the specified coord 

327 # I will check adjacent rings which may include the specified coord 

328 for r in [ringNum - 1, ringNum, ringNum + 1]: 

329 if r < 0 or r >= self.config.numRings: 

330 # Poles will be checked explicitly outside this loop 

331 continue 

332 tractNum = self._raToTractNum(coord.getLongitude(), r) 

333 # Adjacent tracts will also be checked. 

334 for t in [tractNum - 1, tractNum, tractNum + 1]: 

335 # Wrap over raStart 

336 if t < 0: 

337 t = t + self._ringNums[r] 

338 elif t > self._ringNums[r] - 1: 

339 t = t - self._ringNums[r] 

340 

341 extra = 0 

342 if self._version == 0 and t == 0 and r != 0: 

343 # Account for off-by-one error in getRingIndices 

344 # Note that this means that tract 1 gets duplicated. 

345 extra = 1 

346 

347 index = sum(self._ringNums[:r + extra], t + 1) # Allow 1 for south pole 

348 tract = self[index] 

349 if tract.contains(coord): 

350 tractList.append(tract) 

351 

352 # Always check tracts at poles 

353 # Southern cap is 0, Northern cap is the last entry in self 

354 for entry in [0, len(self)-1]: 

355 tract = self[entry] 

356 if tract.contains(coord): 

357 tractList.append(tract) 

358 

359 return tractList 

360 

361 def findTractPatchList(self, coordList): 

362 retList = [] 

363 for coord in coordList: 

364 for tractInfo in self.findAllTracts(coord): 

365 patchList = tractInfo.findPatchList(coordList) 

366 if patchList and not (tractInfo, patchList) in retList: 

367 retList.append((tractInfo, patchList)) 

368 return retList 

369 

370 def updateSha1(self, sha1): 

371 """Add subclass-specific state or configuration options to the SHA1.""" 

372 sha1.update(struct.pack("<id", self.config.numRings, self.config.raStart))