Coverage for python/lsst/pipe/tasks/associationUtils.py: 16%

1# This file is part of pipe_tasks. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22"""Utilities for interfacing with hpgeom. Originally implemented in 

23http://github.com/LSSTDESC/dia_pipe and then translated to hpgeom. 

24""" 

25 

26__all__ = ["toIndex", "toRaDec", "eq2xyz", "eq2xyzVec", "convert_spherical", 

27 "convert_spherical_array", "query_disc", "obj_id_to_ss_object_id", "ss_object_id_to_obj_id"] 

28 

29import hpgeom as hpg 

30import numpy as np 

31 

32 

33def toIndex(nside, ra, dec): 

34 """Return healpix index given ra, dec in degrees 

35 

36 Parameters 

37 ---------- 

38 nside : `int` 

39 Power of 2 nside healpix resolution. 

40 ra : `float` 

41 RA in degrees. 

42 dec : `float` 

43 Declination in degrees 

44 

45 Returns 

46 ------- 

47 index : `int` 

48 Unique healpix pixel ID containing point RA, DEC at resolution nside. 

49 """ 

50 return hpg.angle_to_pixel(nside, ra, dec, nest=False) 

51 

52 

53def toRaDec(nside, index): 

54 """Convert from healpix index to ra,dec in degrees 

55 

56 Parameters 

57 ---------- 

58 nside : `int` 

59 Resolution of healpixel "grid". 

60 index : `int` 

61 Index of the healpix pixel we want to find the location of. 

62 

63 Returns 

64 ------- 

65 pos : `numpy.ndarray`, (2,) 

66 RA and DEC of healpix pixel location in degrees. 

67 """ 

68 ra, dec = hpg.pixel_to_angle(nside, index, nest=False) 

69 return np.dstack((ra, dec))[0] 

70 

71 

72def eq2xyz(ra, dec): 

73 """Convert from equatorial ra,dec in degrees to x,y,z on unit sphere. 

74 

75 Parameters 

76 ---------- 

77 ra : `float` 

78 RA in degrees. 

79 dec : `float` 

80 Declination in degrees 

81 

82 Returns 

83 ------- 

84 xyz : `numpy.ndarray`, (3,) 

85 Float xyz positions on the unit sphere. 

86 """ 

87 phi = np.deg2rad(ra) 

88 theta = np.pi/2 - np.deg2rad(dec) 

89 sintheta = np.sin(theta) 

90 x = sintheta*np.cos(phi) 

91 y = sintheta*np.sin(phi) 

92 z = np.cos(theta) 

93 return np.array([x, y, z]) 

94 

95 

96def eq2xyzVec(ra, dec): 

97 """Convert equatorial ra,dec in degrees to x,y,z on the unit sphere 

98 parameters 

99 

100 Vectorized version of ``eq2xyz`` 

101 

102 Parameters 

103 ---------- 

104 ra : array_like, (N,) 

105 Array of RA in degrees. 

106 dec : array_like, (N,) 

107 Declination in degrees 

108 

109 Returns 

110 ------- 

111 vec : `numpy.ndarray`, (N,3) 

112 Array of unitsphere 3-vectors. 

113 """ 

114 ra = np.array(ra, dtype='f8', ndmin=1, copy=False) 

115 dec = np.array(dec, dtype='f8', ndmin=1, copy=False) 

116 if ra.size != dec.size: 

117 raise ValueError("ra,dec not same size: %s,%s" % (ra.size, dec.size)) 

118 

119 vec = eq2xyz(ra, dec) 

120 

121 return vec 

122 

123 

124def convert_spherical(ra, dec): 

125 """Convert from ra,dec to spherical coordinates. 

126 

127 Used in query_disc. 

128 

129 Parameters 

130 ---------- 

131 ra : `float` 

132 RA in radians. 

133 dec : `float` 

134 Declination in radians 

135 """ 

136 return np.dstack([np.cos(dec*np.pi/180)*np.cos(ra*np.pi/180), 

137 np.cos(dec*np.pi/180)*np.sin(ra*np.pi/180), 

138 np.sin(dec*np.pi/180)])[0] 

139 

140 

141def convert_spherical_array(array): 

142 """Convert from and a array ra,dec to spherical coordinates. 

143 

144 Used in query_disc 

145 

146 Parameters 

147 ---------- 

148 array : `numpy.ndarray`, (N, 2) 

149 (N, 2) Array of RA, DEC values. 

150 

151 Returns 

152 ------- 

153 vecs : `numpy.ndarray`, (N, 3) 

154 Vectors on the unit sphere 

155 """ 

156 ra = array[:, 0] 

157 dec = array[:, 1] 

158 return convert_spherical(ra, dec) 

159 

160 

161def query_disc(nside, ra, dec, max_rad, min_rad=0): 

162 """Get the list of healpix indices within max_rad, min_rad given in radians 

163 around ra,dec given in degrees 

164 

165 Parameters 

166 ---------- 

167 nside : `int` 

168 Resolution of the healpixels to search/return. 

169 ra : `float` 

170 RA in degrees. 

171 dec : `float` 

172 Declination in degrees 

173 max_rad : `float` 

174 Max distance in radians to search nearby healpixels. 

175 min_rad : `float`, optional 

176 Minimum distance in radians to search healpixels. Default = 0. 

177 """ 

178 ra = np.atleast_1d(ra) 

179 dec = np.atleast_1d(dec) 

180 

181 max_rad_deg = np.rad2deg(max_rad) 

182 

183 pixels = np.unique( 

184 [hpg.query_circle(nside, a, b, max_rad_deg, nest=False) 

185 for (a, b) in zip(ra, dec)]) 

186 

187 if min_rad > 0 and len(pixels) > 0: 

188 vec0 = convert_spherical(ra, dec) 

189 min_rad2 = min_rad**2 

190 vecs = convert_spherical_array(toRaDec(nside, pixels)) 

191 dsq = np.sum((vecs - vec0)**2, axis=1) 

192 match = dsq > min_rad2 

193 pixels = pixels[match] 

194 

195 return pixels 

196 

197 

198def obj_id_to_ss_object_id(objID): 

199 """Convert from Minor Planet Center packed provisional object ID to 

200 Rubin ssObjectID. 

201 

202 Parameters 

203 ---------- 

204 objID : `str` 

205 Minor Planet Center packed provisional designation for a small solar 

206 system object. Must be fewer than eight characters. 

207 

208 Returns 

209 ------- 

210 ssObjectID : `int` 

211 Rubin ssObjectID 

212 

213 Raises 

214 ------ 

215 ValueError 

216 Raised if either objID is shorter than 7 or longer than 8 characters or contains 

217 illegal objID characters 

218 """ 

219 if len(objID) > 8: 

220 raise ValueError(f'objID longer than 8 characters: "{objID}"') 

221 if len(objID) < 7: 

222 raise ValueError(f'objID shorter than 7 characters: "{objID}"') 

223 if any([ord(c) > 255 for c in objID]): 

224 raise ValueError(f'{[c for c in objID if ord(c) > 255]} not legal objID characters (ascii [1, 255])') 

225 

226 ssObjectID = ord(objID[0]) 

227 for character in objID[1:]: 

228 ssObjectID <<= 8 

229 ssObjectID += ord(character) 

230 return ssObjectID 

231 

232 

233def ss_object_id_to_obj_id(ssObjectID): 

234 """Convert from Rubin ssObjectID to Minor Planet Center packed provisional 

235 object ID. 

236 

237 Parameters 

238 ---------- 

239 ssObjectID : `int` 

240 Rubin ssObjectID 

241 

242 Returns 

243 ------- 

244 objID : `str` 

245 Minor Planet Center packed provisional designation. 

246 

247 Raises 

248 ------ 

249 """ 

250 if ssObjectID < 0 or ssObjectID >= (1 << 64): 

251 raise ValueError(f'ssObjectID ({ssObjectID}) outside [0, 2^64 - 1].') 

252 

253 objID = ''.join([chr((ssObjectID >> (8 * i)) % 256) for i in reversed(range(0, 8))]) 

254 objID = objID.replace('\x00', '') 

255 return objID