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, flags=0): 

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 flags : `int`, optional 

208 Eight free bits to enable future decoupling between Minor Planet Center 

209 and Rubin. Zero by default, should not be changed unless we need to 

210 move away from a 1:1 mapping with the MPC. Must be within [0, 255]. 

211 

212 Returns 

213 ------- 

214 ssObjectID : `int` 

215 Rubin ssObjectID 

216 

217 Raises 

218 ------ 

219 ValueError 

220 Raised if either objID is longer than 7 characters or flags is greater 

221 than 255 or less than 0. 

222 """ 

223 if len(objID) > 7: 

224 raise ValueError(f'objID longer than 7 characters: "{objID}"') 

225 if len(objID) < 7: 

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

227 if flags < 0 or flags > 255: 

228 raise ValueError(f'Flags ({flags}) outside [0, 255].') 

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

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

231 

232 ssObjectID = flags 

233 for character in objID: 

234 ssObjectID <<= 8 

235 ssObjectID += ord(character) 

236 return ssObjectID 

237 

238 

239def ss_object_id_to_obj_id(ssObjectID): 

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

241 object ID. 

242 

243 Parameters 

244 ---------- 

245 ssObjectID : `int` 

246 Rubin ssObjectID 

247 

248 Returns 

249 ------- 

250 objID : `str` 

251 Minor Planet Center packed provisional designation. 

252 

253 flags : `int` 

254 Rubin flags (not yet defined, but usable in case we decouple from MPC). 

255 

256 Raises 

257 ------ 

258 """ 

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

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

261 

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

263 return objID, ssObjectID >> (8 * 7) % 256