1 """module for coordinate manipulation (conversions, calculations etc.)
2
3 (c) 2007-2012 Matt Hilton
4
5 (c) 2013-2014 Matt Hilton & Steven Boada
6
7 U{http://astlib.sourceforge.net}
8
9 """
10
11 import numpy
12 from PyWCSTools import wcscon
13
14
15
17 """Converts a delimited string of Hours:Minutes:Seconds format into decimal
18 degrees.
19
20 @type RAString: string
21 @param RAString: coordinate string in H:M:S format
22 @type delimiter: string
23 @param delimiter: delimiter character in RAString
24 @rtype: float
25 @return: coordinate in decimal degrees
26
27 """
28
29 if delimiter == "":
30 RABits = str(RAString).split()
31 else:
32 RABits = str(RAString).split(delimiter)
33 if len(RABits) > 1:
34 RAHDecimal = float(RABits[0])
35 if len(RABits) > 1:
36 RAHDecimal = RAHDecimal + (float(RABits[1]) / 60.0)
37 if len(RABits) > 2:
38 RAHDecimal = RAHDecimal + (float(RABits[2]) / 3600.0)
39 RADeg = (RAHDecimal / 24.0) * 360.0
40 else:
41 RADeg = float(RAString)
42
43 return RADeg
44
45
46
48 """Converts a delimited string of Degrees:Minutes:Seconds format into
49 decimal degrees.
50
51 @type decString: string
52 @param decString: coordinate string in D:M:S format
53 @type delimiter: string
54 @param delimiter: delimiter character in decString
55 @rtype: float
56 @return: coordinate in decimal degrees
57
58 """
59
60 if delimiter == "":
61 decBits = str(decString).split()
62 else:
63 decBits = str(decString).split(delimiter)
64 if len(decBits) > 1:
65 decDeg = float(decBits[0])
66 if decBits[0].find("-") != -1:
67 if len(decBits) > 1:
68 decDeg = decDeg - (float(decBits[1]) / 60.0)
69 if len(decBits) > 2:
70 decDeg = decDeg - (float(decBits[2]) / 3600.0)
71 else:
72 if len(decBits) > 1:
73 decDeg = decDeg + (float(decBits[1]) / 60.0)
74 if len(decBits) > 2:
75 decDeg = decDeg + (float(decBits[2]) / 3600.0)
76 else:
77 decDeg = float(decString)
78
79 return decDeg
80
81
82
84 """Converts decimal degrees to string in Hours:Minutes:Seconds format with
85 user specified delimiter.
86
87 @type RADeg: float
88 @param RADeg: coordinate in decimal degrees
89 @type delimiter: string
90 @param delimiter: delimiter character in returned string
91 @rtype: string
92 @return: coordinate string in H:M:S format
93
94 """
95 hours = (RADeg / 360.0) * 24
96
97 if 1 <= hours < 10:
98 sHours = "0" + str(hours)[0]
99 elif hours >= 10:
100 sHours = str(hours)[:2]
101 elif hours < 1:
102 sHours = "00"
103
104 if str(hours).find(".") == -1:
105 mins = float(hours) * 60.0
106 else:
107 mins = float(str(hours)[str(hours).index("."):]) * 60.0
108
109 if 1 <= mins < 10:
110 sMins = "0" + str(mins)[:1]
111 elif mins >= 10:
112 sMins = str(mins)[:2]
113 elif mins < 1:
114 sMins = "00"
115
116 secs = (hours - (float(sHours) + float(sMins) / 60.0)) * 3600.0
117
118 if 0.001 < secs < 10:
119 sSecs = "0" + str(secs)[:str(secs).find(".") + 4]
120 elif secs < 0.0001:
121 sSecs = "00.001"
122 else:
123 sSecs = str(secs)[:str(secs).find(".") + 4]
124 if len(sSecs) < 5:
125 sSecs = sSecs + "00"
126
127 if float(sSecs) == 60.000:
128 sSecs = "00.00"
129 sMins = str(int(sMins) + 1)
130 if int(sMins) == 60:
131 sMins = "00"
132
133
134 return sHours + delimiter + sMins + delimiter + sSecs
135
136
137
139 """Converts decimal degrees to string in Degrees:Minutes:Seconds format
140 with user specified delimiter.
141
142 @type decDeg: float
143 @param decDeg: coordinate in decimal degrees
144 @type delimiter: string
145 @param delimiter: delimiter character in returned string
146 @rtype: string
147 @return: coordinate string in D:M:S format
148
149 """
150
151 if decDeg > 0:
152
153 if 1 <= decDeg < 10:
154 sDeg = "0" + str(decDeg)[0]
155 elif decDeg >= 10:
156 sDeg = str(decDeg)[:2]
157 elif decDeg < 1:
158 sDeg = "00"
159
160 if str(decDeg).find(".") == -1:
161 mins = float(decDeg) * 60.0
162 else:
163 mins = float(str(decDeg)[str(decDeg).index("."):]) * 60
164
165 if 1 <= mins < 10:
166 sMins = "0" + str(mins)[:1]
167 elif mins >= 10:
168 sMins = str(mins)[:2]
169 elif mins < 1:
170 sMins = "00"
171
172 secs = (decDeg - (float(sDeg) + float(sMins) / 60.0)) * 3600.0
173
174 if 0 < secs < 10:
175 sSecs = "0" + str(secs)[:str(secs).find(".") + 3]
176 elif secs < 0.001:
177 sSecs = "00.00"
178 else:
179 sSecs = str(secs)[:str(secs).find(".") + 3]
180 if len(sSecs) < 5:
181 sSecs = sSecs + "0"
182
183 if float(sSecs) == 60.00:
184 sSecs = "00.00"
185 sMins = str(int(sMins) + 1)
186 if int(sMins) == 60:
187 sMins = "00"
188 sDeg = str(int(sDeg) + 1)
189
190 return "+" + sDeg + delimiter + sMins + delimiter + sSecs
191
192 else:
193
194 if -10 < decDeg <= -1:
195 sDeg = "-0" + str(decDeg)[1]
196 elif decDeg <= -10:
197 sDeg = str(decDeg)[:3]
198 elif decDeg > -1:
199 sDeg = "-00"
200
201 if str(decDeg).find(".") == -1:
202 mins = float(decDeg) * -60.0
203 else:
204 mins = float(str(decDeg)[str(decDeg).index("."):]) * 60
205
206 if 1 <= mins < 10:
207 sMins = "0" + str(mins)[:1]
208 elif mins >= 10:
209 sMins = str(mins)[:2]
210 elif mins < 1:
211 sMins = "00"
212
213 secs = (decDeg - (float(sDeg) - float(sMins) / 60.0)) * 3600.0
214
215
216 if -10 < secs < 0:
217 sSecs = "0" + str(secs)[1:str(secs).find(".") + 3]
218 elif secs > -0.001:
219 sSecs = "00.00"
220 else:
221 sSecs = str(secs)[1:str(secs).find(".") + 3]
222 if len(sSecs) < 5:
223 sSecs = sSecs + "0"
224
225 if float(sSecs) == 60.00:
226 sSecs = "00.00"
227 sMins = str(int(sMins) + 1)
228 if int(sMins) == 60:
229 sMins = "00"
230 sDeg = str(int(sDeg) - 1)
231
232 return sDeg + delimiter + sMins + delimiter + sSecs
233
234
235
237 """
238 Convert Equatorial coordinates to Cartesian coordinates. Return a tuple
239 (x, y, z) in the same unit of the input distance. This is the inverse of
240 cart2eq.
241
242 @type RA: float
243 @param RA: Right Ascension in decimal degrees
244 @type DEC: float
245 @param DEC: Declination in decimal degrees
246 @type r: float
247 @param r: Distance to the object.
248 @rtype: tuple
249 @return: Tuple of (x, y, z) in same unit as the input distance.
250
251 """
252
253 RA = numpy.radians(RA)
254 DEC = numpy.radians(DEC)
255
256 x = r * numpy.cos(RA) * numpy.cos(DEC)
257 y = r * numpy.sin(RA) * numpy.cos(DEC)
258 z = r * numpy.sin(DEC)
259
260 return x, y, z
261
262
263
265 """
266 Convert Cartesian coordinates to Equatorial coordinates. Returns a tuple of
267 (RA, DEC, r), with RA and DEC given in decimal degrees and r in the same
268 unit as the input.
269
270 @type x: float
271 @param x: x coordinate
272 @type y: float
273 @param y: x coordinate
274 @type z: float
275 @param z: x coordinate
276 @rtype: tuple
277 @return: Tuple of (RA, DEC, r)
278
279 """
280
281 r = numpy.sqrt(numpy.power(x, 2) + numpy.power(y, 2) + numpy.power(z, 2))
282 ra = numpy.arctan2(y, x)
283 try:
284 ra[ra < 0] = 2.0 * numpy.pi + ra[ra < 0]
285 except TypeError:
286 ra = ra if ra >= 0 else (2.0 * numpy.pi + ra)
287 dec = numpy.arcsin(z / r)
288
289 ra = numpy.degrees(ra)
290 dec = numpy.degrees(dec)
291
292 return ra, dec, r
293
294
295
297 """Calculates the angular separation of two positions on the sky (specified
298 in decimal degrees) in decimal degrees, assuming a tangent plane projection
299 (so separation has to be <90 deg). Note that RADeg2, decDeg2 can be numpy
300 arrays.
301
302 @type RADeg1: float
303 @param RADeg1: R.A. in decimal degrees for position 1
304 @type decDeg1: float
305 @param decDeg1: dec. in decimal degrees for position 1
306 @type RADeg2: float or numpy array
307 @param RADeg2: R.A. in decimal degrees for position 2
308 @type decDeg2: float or numpy array
309 @param decDeg2: dec. in decimal degrees for position 2
310 @rtype: float or numpy array, depending upon type of RADeg2, decDeg2
311 @return: angular separation in decimal degrees
312
313 """
314 cRA = numpy.radians(RADeg1)
315 cDec = numpy.radians(decDeg1)
316
317 gRA = numpy.radians(RADeg2)
318 gDec = numpy.radians(decDeg2)
319
320
321
322 cosC = ((numpy.sin(gDec) * numpy.sin(cDec)) +
323 (numpy.cos(gDec) * numpy.cos(cDec) * numpy.cos(gRA - cRA)))
324 x = (numpy.cos(cDec) * numpy.sin(gRA - cRA)) / cosC
325 y = (((numpy.cos(gDec) * numpy.sin(cDec)) -
326 (numpy.sin(gDec) * numpy.cos(cDec) * numpy.cos(gRA - cRA))) / cosC)
327 r = numpy.degrees(numpy.sqrt(x * x + y * y))
328
329 return r
330
331
332
334 """Computes new right ascension and declination shifted from the original
335 by some delta RA and delta DEC. Input position is decimal degrees. Shifts
336 (deltaRA, deltaDec) are arcseconds, and output is decimal degrees. Based on
337 an IDL routine of the same name.
338
339 @param ra1: float
340 @type ra1: R.A. in decimal degrees
341 @param dec1: float
342 @type dec1: dec. in decimal degrees
343 @param deltaRA: float
344 @type deltaRA: shift in R.A. in arcseconds
345 @param deltaDec: float
346 @type deltaDec: shift in dec. in arcseconds
347 @rtype: float [newRA, newDec]
348 @return: shifted R.A. and dec.
349
350 """
351
352 d2r = numpy.pi / 180.
353 as2r = numpy.pi / 648000.
354
355
356
357 dcrad1 = dec1 * d2r
358 shiftRArad = deltaRA * as2r
359
360
361
362
363 sindis = numpy.sin(shiftRArad / 2.0)
364 sindelRA = sindis / numpy.cos(dcrad1)
365 delra = 2.0 * numpy.arcsin(sindelRA) / d2r
366
367
368 ra2 = ra1 + delra
369 dec2 = dec1 + deltaDec / 3600.0
370
371
372 if ra2 > 360:
373 ra2 -= 360
374 elif ra2 < 0:
375 ra2 += 360
376
377 return ra2, dec2
378
379
380
381 -def convertCoords(inputSystem, outputSystem, coordX, coordY, epoch):
382 """Converts specified coordinates (given in decimal degrees) between J2000,
383 B1950, and Galactic.
384
385 @type inputSystem: string
386 @param inputSystem: system of the input coordinates (either "J2000",
387 "B1950" or "GALACTIC")
388 @type outputSystem: string
389 @param outputSystem: system of the returned coordinates (either "J2000",
390 "B1950" or "GALACTIC")
391 @type coordX: float
392 @param coordX: longitude coordinate in decimal degrees, e.g. R. A.
393 @type coordY: float
394 @param coordY: latitude coordinate in decimal degrees, e.g. dec.
395 @type epoch: float
396 @param epoch: epoch of the input coordinates
397 @rtype: list
398 @return: coordinates in decimal degrees in requested output system
399
400 """
401
402 if inputSystem == "J2000" or inputSystem == "B1950" or inputSystem == "GALACTIC":
403 if outputSystem == "J2000" or outputSystem == "B1950" or \
404 outputSystem == "GALACTIC":
405
406 outCoords = wcscon.wcscon(
407 wcscon.wcscsys(inputSystem), wcscon.wcscsys(outputSystem), 0,
408 0, coordX, coordY, epoch)
409
410 return outCoords
411
412 raise Exception("inputSystem and outputSystem must be 'J2000', 'B1950'"
413 "or 'GALACTIC'")
414
415
416
418 """ Calculates the area of the quadrangle on the sky given by the
419 intersection of a right ascension lune and declination zone. By default,
420 the area is returned in square degrees.
421
422 @type RA1: float
423 @param RA1: Right ascension in decimal degrees
424 @type RA2: float
425 @param RA2: Right ascension in decimal degrees
426 @type DEC1: float
427 @param DEC1: Declination in decimal degrees
428 @type DEC2: float
429 @param DEC2: Declination in decimal degrees
430 @type units: bool
431 @param units: True: Area is given in square degrees. False: Area is given
432 as a fraction of the total sky.
433
434 """
435
436
437 RA1 = numpy.radians(RA1)
438 RA2 = numpy.radians(RA2)
439 DEC1 = numpy.radians(DEC1)
440 DEC2 = numpy.radians(DEC2)
441
442 if units:
443 return abs(RA1 - RA2) * abs(numpy.sin(DEC1) - numpy.sin(DEC2)) *\
444 (180 / numpy.pi)**2
445 else:
446 return abs(RA1 - RA2) * abs(numpy.sin(DEC1) - numpy.sin(DEC2))
447
448
449
451 """Calculates minimum and maximum RA, dec coords needed to define a box
452 enclosing a circle of radius radiusSkyDeg around the given RADeg, decDeg
453 coordinates. Useful for freeform queries of e.g. SDSS, UKIDSS etc.. Uses
454 L{calcAngSepDeg}, so has the same limitations.
455
456 @type RADeg: float
457 @param RADeg: RA coordinate of centre of search region
458 @type decDeg: float
459 @param decDeg: dec coordinate of centre of search region
460 @type radiusSkyDeg: float
461 @param radiusSkyDeg: radius in degrees on the sky used to define search
462 region
463 @rtype: list
464 @return: [RAMin, RAMax, decMin, decMax] - coordinates in decimal degrees
465 defining search box
466
467 """
468
469 tolerance = 1e-5
470 targetHalfSizeSkyDeg = radiusSkyDeg
471 funcCalls = ["calcAngSepDeg(RADeg, decDeg, guess, decDeg)",
472 "calcAngSepDeg(RADeg, decDeg, guess, decDeg)",
473 "calcAngSepDeg(RADeg, decDeg, RADeg, guess)",
474 "calcAngSepDeg(RADeg, decDeg, RADeg, guess)"]
475 coords = [RADeg, RADeg, decDeg, decDeg]
476 signs = [1.0, -1.0, 1.0, -1.0]
477 results = []
478 for f, c, sign in zip(funcCalls, coords, signs):
479
480 maxGuess = sign * targetHalfSizeSkyDeg * 10.0
481 minGuess = sign * targetHalfSizeSkyDeg / 10.0
482 guessStep = (maxGuess - minGuess) / 10.0
483 guesses = numpy.arange(minGuess + c, maxGuess + c, guessStep)
484 for i in range(50):
485 minSizeDiff = 1e6
486 bestGuess = None
487 for guess in guesses:
488 sizeDiff = abs(eval(f) - targetHalfSizeSkyDeg)
489 if sizeDiff < minSizeDiff:
490 minSizeDiff = sizeDiff
491 bestGuess = guess
492 if minSizeDiff < tolerance:
493 break
494 else:
495 guessRange = abs((maxGuess - minGuess))
496 maxGuess = bestGuess + guessRange / 4.0
497 minGuess = bestGuess - guessRange / 4.0
498 guessStep = (maxGuess - minGuess) / 10.0
499 guesses = numpy.arange(minGuess, maxGuess, guessStep)
500 results.append(bestGuess)
501
502 RAMax = results[0]
503 RAMin = results[1]
504 decMax = results[2]
505 decMin = results[3]
506
507 return [RAMin, RAMax, decMin, decMax]
508
509
511 """
512 Make Aitoff map projection.
513
514 Take traditional longitude and latitude in radians and return a
515 tuple (x, y).
516
517 Notice that traditionally longitude is in [-pi:pi] from the meridian,
518 and latitude is in [-pi/2:pi/2] from the equator. So, for example, if
519 you would like to make a galactic map projection centered on the galactic
520 center, before passing galactic longitude l to the function you should
521 first do:
522 l = l if l <= numpy.pi else l - 2 * numpy.pi
523
524 Keyword arguments:
525 lon -- Traditional longitude in radians, in range [-pi:pi]
526 lat -- Traditional latitude in radians, in range [-pi/2:pi/2]
527 """
528
529 def sinc(x):
530
531 if not x:
532 return 0
533 else:
534 return numpy.sin(x) / x
535
536 x = numpy.zeros_like(lon)
537 y = numpy.zeros_like(lat)
538
539
540 if lon > numpy.pi or lon < -numpy.pi or lat > numpy.pi / 2 or \
541 lat < -numpy.pi / 2:
542 print('Aitoff: Input longitude and latitude out of range.\n')
543 print(' lon: [-pi,pi]; lat: [-pi/2,pi/2].\n')
544 return None
545
546
547
548 if lon == 0 and lat == 0:
549 return 0.0, 0.0
550
551 alpha = numpy.acos(numpy.cos(lat) * numpy.cos(lon / 2.0))
552
553
554 x = 2.0 * numpy.cos(lat) * numpy.sin(lon / 2.0) / sinc(alpha)
555
556 y = numpy.sin(lat) / sinc(alpha)
557
558 return x, y
559