1#!/usr/local/bin/python 

2# encoding: utf-8 

3""" 

4*Tools for working with Hierarchical Triangular Meshes, including coordinate crossmatching* 

5 

6:Author: 

7 David Young (originally forked from Erin Sheldon's esutil - esheldon) 

8""" 

9from __future__ import absolute_import 

10from __future__ import division 

11from past.utils import old_div 

12import sys 

13import os 

14os.environ['TERM'] = 'vt100' 

15from fundamentals import tools 

16from astrocalc.coords import unit_conversion 

17from . import _htmcCode 

18import numpy 

19from sys import stdout 

20 

21class HTM(_htmcCode.HTMC): 

22 """ 

23 *A Hierarchical Triangular Mesh object* 

24 

25 **Key Arguments** 

26 

27 - ``depth`` -- the depth of the mesh you wish to create. Default *16* 

28  

29 

30 **Usage** 

31 

32 To generate a mesh object: 

33 

34 ```python 

35 from HMpTy import HTM 

36 mesh16 = HTM( 

37 depth=16 

38 ) 

39 ``` 

40  

41 """ 

42 

43 @property 

44 def depth( 

45 self): 

46 """*the depth of the HTM tree* 

47 

48 **Usage** 

49 

50 ```python 

51 mesh.depth 

52 ``` 

53  

54 """ 

55 return super(HTM, self).depth() 

56 

57 @property 

58 def area( 

59 self): 

60 """*The mean area of triangles in this mesh in units of square degrees.* 

61 

62 **Usage** 

63 

64 ```python 

65 mesh.area 

66 ``` 

67  

68 """ 

69 pi = numpy.pi 

70 area0 = 4.0 * pi / 8.0 

71 areadiv = 4.0 ** self.depth 

72 area = old_div(area0, areadiv) * (180.0 / pi) ** 2 

73 return area 

74 

75 def lookup_id( 

76 self, 

77 ra, 

78 dec): 

79 """*Lookup the ID of HTM trixel that a coordinate or lists of coordinates lie on* 

80 

81 **Key Arguments** 

82 

83 - ``ra`` -- list, numpy array or single ra value (first coordinate set) 

84 - ``dec`` -- list, numpy array or single dec value (first coordinate set - must match ra1 array length) 

85  

86 

87 **Return** 

88 

89 - ``htmIds`` -- a list of HTM trixel ids the coordinates lie on 

90  

91 

92 **Usage** 

93 

94 To find the trixel IDs that a set of coordinate lie on: 

95 

96 ```python 

97 raList1 = ["13:20:00.00", 200.0, "13:20:00.00", 175.23, 21.36] 

98 decList1 = ["+24:18:00.00", 24.3, "+24:18:00.00", -28.25, -15.32] 

99 

100 htmids = mesh.lookup_id(raList1, decList1) 

101 for h, r, d in zip(htmids, raList1, decList1): 

102 print(r, d, " --> ", h) 

103 ``` 

104  

105 """ 

106 self.log.debug('starting the ``lookup_id`` method') 

107 

108 from astrocalc.coords import coordinates_to_array 

109 raArray, decArray = coordinates_to_array( 

110 log=self.log, 

111 ra=ra, 

112 dec=dec 

113 ) 

114 

115 self.log.debug('completed the ``lookup_id`` method') 

116 return super(HTM, self).lookup_id(raArray, decArray) 

117 

118 # use the tab-trigger below for new method 

119 # xt-class-method 

120 

121 def intersect(self, ra, dec, radius, inclusive=True, convertCoordinates=True): 

122 """*return IDs of all triangles contained within and/or intersecting a circle centered on a given ra and dec* 

123 

124 **Key Arguments** 

125 

126 - ``ra`` -- RA of central point in decimal degrees or sexagesimal 

127 - ``dec`` -- DEC of central point in decimal degrees or sexagesimal 

128 - ``radius`` -- radius of circle in degrees 

129 - ``inclusive`` -- include IDs of triangles that intersect the circle as well as those completely inclosed by the circle. Default *True* 

130 - ``convertCoordinates`` -- convert the corrdinates passed to intersect. Default *True* 

131 - 

132  

133 

134 **Return** 

135 

136 - ``trixelArray`` -- a numpy array of the match trixel IDs 

137  

138 

139 **Usage** 

140 

141 To return the trixels overlapping a circle with a 10 arcsec radius centred at 23:25:53.56, +26:54:23.9 

142 

143 ```python 

144 overlappingTrixels = mesh16.intersect( 

145 ra="23:25:53.56", 

146 dec="+26:54:23.9", 

147 radius=10 / (60 * 60), 

148 inclusive=True 

149 ) 

150 ``` 

151 

152 Or to return the trixels completing enclosed by a circle with a 1 degree radius centred at 23:25:53.56, +26:54:23.9 

153 

154 ```python 

155 overlappingTrixels = mesh16.intersect( 

156 ra="23:25:53.56", 

157 dec="+26:54:23.9", 

158 radius=1, 

159 inclusive=False 

160 ) 

161 ``` 

162  

163 """ 

164 # CONVERT RA AND DEC DECIMAL DEGREES 

165 

166 if convertCoordinates == True: 

167 converter = unit_conversion( 

168 log=self.log 

169 ) 

170 ra = converter.ra_sexegesimal_to_decimal( 

171 ra=ra 

172 ) 

173 dec = converter.dec_sexegesimal_to_decimal( 

174 dec=dec 

175 ) 

176 

177 if inclusive: 

178 inc = 1 

179 else: 

180 inc = 0 

181 return super(HTM, self).intersect(ra, dec, radius, inc) 

182 

183 def match(self, ra1, dec1, ra2, dec2, radius, maxmatch=1, convertToArray=True): 

184 """*Crossmatch two lists of ra/dec points* 

185 

186 This is very efficient for large search angles and large lists. Note, if you need to match against the same points many times, you should use a `Matcher` object 

187 

188 **Key Arguments** 

189 

190 - ``ra1`` -- list, numpy array or single ra value (first coordinate set) 

191 - ``dec1`` -- list, numpy array or single dec value (first coordinate set - must match ra1 array length) 

192 - ``ra2`` -- list, numpy array or single ra value (second coordinate set) 

193 - ``dec2`` -- list, numpy array or single dec value (second coordinate set - must match ra2 array length) 

194 - ``radius`` -- search radius in degrees. Can be list, numpy array or single value. If list or numpy array must be same length as ra1 array length) 

195 - ``maxmatch`` -- maximum number of matches to return. Set to `0` to match all points. Default *1* (i.e. closest match) 

196 - ``convertToArray`` -- convert the coordinates into an array. Default *True*. Can bypass the conversion check if you are sure coordinates in numpy array 

197  

198 

199 **Return** 

200 

201 - ``matchIndices1`` -- match indices for list1 (ra1, dec1) 

202 - ``matchIndices2`` -- match indices for list2 (ra2, dec2) 

203 - ``sepDeg`` -- separations between matched corrdinates in degrees. All returned arrays are the same size 

204  

205 

206 **Usage** 

207 

208 To match 2 lists of corrdinates try something like this: 

209 

210 ```python 

211 twoArcsec = 2.0 / 3600. 

212 raList1 = [200.0, 200.0, 200.0, 175.23, 21.36] 

213 decList1 = [24.3, 24.3, 24.3, -28.25, -15.32] 

214 raList2 = [200.0, 200.0, 200.0, 175.23, 55.25] 

215 decList2 = [24.3 + 0.75 * twoArcsec, 24.3 + 0.25 * twoArcsec, 

216 24.3 - 0.33 * twoArcsec, -28.25 + 0.58 * twoArcsec, 75.22] 

217 matchIndices1, matchIndices2, seps = mesh.match( 

218 ra1=raList1, 

219 dec1=decList1, 

220 ra2=raList2, 

221 dec2=decList2, 

222 radius=twoArcsec, 

223 maxmatch=0 

224 ) 

225 

226 for m1, m2, s in zip(matchIndices1, matchIndices2, seps): 

227 print(raList1[m1], decList1[m1], " -> ", s * 3600., " arcsec -> ", raList2[m2], decList2[m2]) 

228 ``` 

229 

230 Note from the print statement, you can index the arrays ``raList1``, ``decList1`` with the ``matchIndices1`` array values and ``raList2``, ``decList2`` with the ``matchIndices2`` values. 

231  

232 """ 

233 

234 # CONVERT LISTS AND SINGLE VALUES TO ARRAYS OF FLOATS 

235 ra1 = numpy.array(ra1, dtype='f8', ndmin=1, copy=False) 

236 dec1 = numpy.array(dec1, dtype='f8', ndmin=1, copy=False) 

237 ra2 = numpy.array(ra2, dtype='f8', ndmin=1, copy=False) 

238 dec2 = numpy.array(dec2, dtype='f8', ndmin=1, copy=False) 

239 radius = numpy.array(radius, dtype='f8', ndmin=1, copy=False) 

240 

241 # CHECK ARRAY SIZES MATCH 

242 if (ra1.size != dec1.size 

243 or ra2.size != ra2.size): 

244 stup = (ra1.size, dec1.size, ra2.size, dec2.size) 

245 raise ValueError("ra1 must equal dec1 in size " 

246 "and ra2 must equal dec2 in size, " 

247 "got %d,%d and %d,%d" % stup) 

248 

249 if radius.size != 1 and radius.size != ra2.size: 

250 raise ValueError("radius size (%d) != 1 and" 

251 " != ra2,dec2 size (%d)" % (radius.size, ra2.size)) 

252 

253 # QUICK TRIMMING IN DEC SPACE OF BOTH SETS OF ARRAYS 

254 decMatchIndices2 = (numpy.abs(dec1[:, None] - dec2) < radius).any(0) 

255 decMatchIndices2 = numpy.where(decMatchIndices2)[0] 

256 ra2a = ra2[decMatchIndices2] 

257 dec2a = dec2[decMatchIndices2] 

258 

259 decMatchIndices1 = (numpy.abs(dec2[:, None] - dec1) < radius).any(0) 

260 decMatchIndices1 = numpy.where(decMatchIndices1)[0] 

261 ra1a = ra1[decMatchIndices1] 

262 dec1a = dec1[decMatchIndices1] 

263 

264 # new way using a Matcher 

265 depth = self.depth 

266 matcher = Matcher( 

267 log=self.log, 

268 depth=depth, 

269 ra=ra2a, 

270 dec=dec2a, 

271 convertToArray=convertToArray) 

272 matchIndices1, matchIndices2, seps = matcher.match( 

273 ra=ra1a, 

274 dec=dec1a, 

275 radius=radius, 

276 maxmatch=maxmatch) 

277 

278 matchIndices1 = decMatchIndices1[matchIndices1] 

279 matchIndices2 = decMatchIndices2[matchIndices2] 

280 

281 return matchIndices1, matchIndices2, seps 

282 

283class Matcher(_htmcCode.Matcher): 

284 """*A matcher-array object to match other arrays of ra,dec against* 

285 

286 The Matcher object is initialized with a set of ra,dec coordinates and can then be used and reused to match against other sets of coordinates 

287 

288 **Key Arguments** 

289 

290 - ``log`` -- logger 

291 - ``depth`` -- the depth of the mesh generate the Matcher object at. Default *16* 

292 - ``ra`` -- list, numpy array or single ra value 

293 - ``dec`` -- --list, numpy array or single dec value (must match ra array length) 

294 - ``convertToArray`` -- convert the coordinates into an array. Default *True*. Can bypass the conversion check if you are sure coordinates in numpy array 

295  

296 

297 **Return** 

298 

299 - None 

300  

301 

302 **Usage** 

303 

304 If we have a set of coordinates such that: 

305 

306 ```python 

307 raList1 = [200.0, 200.0, 200.0, 175.23, 21.36] 

308 decList1 = [24.3, 24.3, 24.3, -28.25, -15.32] 

309 ``` 

310 

311 We can initialise a matcher object like so: 

312 

313 ```python 

314 from HMpTy import Matcher 

315 coordinateSet = Matcher( 

316 log=log, 

317 ra=raList1, 

318 dec=decList1, 

319 depth=16 

320 ) 

321 ``` 

322  

323 """ 

324 

325 def __init__( 

326 self, 

327 ra, 

328 dec, 

329 depth=16, 

330 log=False, 

331 convertToArray=True): 

332 

333 self.convertToArray = convertToArray 

334 

335 if log == False: 

336 if log == False: 

337 from fundamentals.logs import emptyLogger 

338 self.log = emptyLogger() 

339 else: 

340 self.log = log 

341 

342 if convertToArray == True: 

343 from astrocalc.coords import coordinates_to_array 

344 ra, dec = coordinates_to_array( 

345 log=log, 

346 ra=ra, 

347 dec=dec 

348 ) 

349 

350 if ra.size != dec.size: 

351 raise ValueError("ra size (%d) != " 

352 "dec size (%d)" % (ra.size, dec.size)) 

353 

354 super(Matcher, self).__init__(depth, ra, dec) 

355 

356 @property 

357 def depth( 

358 self): 

359 """*the depth of the Matcher object* 

360 

361 **Usage** 

362 

363 ```python 

364 coordinateSet.depth 

365 ``` 

366  

367 """ 

368 return super(Matcher, self).depth 

369 

370 def match(self, ra, dec, radius, maxmatch=1): 

371 """*match a corrdinate set against this Matcher object's coordinate set* 

372 

373 **Key Arguments** 

374 

375 - ``ra`` -- list, numpy array or single ra value 

376 - ``dec`` -- --list, numpy array or single dec value (must match ra array length) 

377 - ``radius`` -- radius of circle in degrees 

378 - ``maxmatch`` -- maximum number of matches to return. Set to `0` to match all points. Default *1* (i.e. closest match) 

379  

380 

381 **Return** 

382 

383 - None 

384  

385 

386 **Usage** 

387 

388 Once we have initialised a Matcher coordinateSet object we can match other coordinate sets against it: 

389 

390 ```python 

391 twoArcsec = 2.0 / 3600. 

392 raList2 = [200.0, 200.0, 200.0, 175.23, 55.25] 

393 decList2 = [24.3 + 0.75 * twoArcsec, 24.3 + 0.25 * twoArcsec, 

394 24.3 - 0.33 * twoArcsec, -28.25 + 0.58 * twoArcsec, 75.22] 

395 

396 matchIndices1, matchIndices2, seps = coordinateSet.match( 

397 ra=raList2, 

398 dec=decList2, 

399 radius=twoArcsec, 

400 maxmatch=0 

401 ) 

402 

403 for m1, m2, s in zip(matchIndices1, matchIndices2, seps): 

404 print(raList1[m1], decList1[m1], " -> ", s * 3600., " arcsec -> ", raList2[m2], decList2[m2]) 

405 ``` 

406 

407 Or to return just the nearest matches: 

408 

409 ```python 

410 matchIndices1, matchIndices2, seps = coordinateSet.match( 

411 ra=raList2, 

412 dec=decList2, 

413 radius=twoArcsec, 

414 maxmatch=1 

415 ) 

416 ``` 

417 

418 Note from the print statement, you can index the arrays ``raList1``, ``decList1`` with the ``matchIndices1`` array values and ``raList2``, ``decList2`` with the ``matchIndices2`` values. 

419  

420 """ 

421 

422 if self.convertToArray == True: 

423 from astrocalc.coords import coordinates_to_array 

424 ra, dec = coordinates_to_array( 

425 log=self.log, 

426 ra=ra, 

427 dec=dec 

428 ) 

429 

430 radius = numpy.array(radius, dtype='f8', ndmin=1, copy=False) 

431 

432 if ra.size != dec.size: 

433 raise ValueError("ra size (%d) != " 

434 "dec size (%d)" % (ra.size, dec.size)) 

435 

436 if radius.size != 1 and radius.size != ra.size: 

437 raise ValueError("radius size (%d) != 1 and" 

438 " != ra,dec size (%d)" % (radius.size, ra.size)) 

439 

440 return super(Matcher, self).match(ra, dec, radius, maxmatch, False)