Project

General

Profile

Download (15.6 KB) Statistics
| Branch: | Revision:

git_sitools_idoc / common / table_creation / astCoords.py @ 97d86a8f

1 97d86a8f Alessandro_N
"""module for coordinate manipulation (conversions, calculations etc.)
2

3
(c) 2007-2012 Matt Hilton
4

5
(c) 2013 Matt Hilton & Steven Boada
6

7
U{http://astlib.sourceforge.net}
8

9
"""
10
11
import sys
12
import math
13
import numpy
14
import string
15
from PyWCSTools import wcscon
16
17
#-----------------------------------------------------------------------------
18
def hms2decimal(RAString, delimiter):
19
    """Converts a delimited string of Hours:Minutes:Seconds format into decimal
20
    degrees.
21

22
    @type RAString: string
23
    @param RAString: coordinate string in H:M:S format
24
    @type delimiter: string
25
    @param delimiter: delimiter character in RAString
26
    @rtype: float
27
    @return: coordinate in decimal degrees
28

29
    """
30
    # is it in HH:MM:SS format?
31
    if delimiter == "":
32
        RABits = str(RAString).split()
33
    else:
34
        RABits = str(RAString).split(delimiter)
35
    if len(RABits) > 1:
36
        RAHDecimal = float(RABits[0])
37
        if len(RABits) > 1:
38
            RAHDecimal = RAHDecimal+(float(RABits[1])/60.0)
39
        if len(RABits) > 2:
40
            RAHDecimal = RAHDecimal+(float(RABits[2])/3600.0)
41
        RADeg = (RAHDecimal/24.0)*360.0
42
    else:
43
        RADeg = float(RAString)
44
45
    return RADeg
46
47
#-----------------------------------------------------------------------------
48
def dms2decimal(decString, delimiter):
49
    """Converts a delimited string of Degrees:Minutes:Seconds format into
50
    decimal degrees.
51

52
    @type decString: string
53
    @param decString: coordinate string in D:M:S format
54
    @type delimiter: string
55
    @param delimiter: delimiter character in decString
56
    @rtype: float
57
    @return: coordinate in decimal degrees
58

59
    """
60
    # is it in DD:MM:SS format?
61
    if delimiter == "":
62
        decBits = str(decString).split()
63
    else:
64
        decBits = str(decString).split(delimiter)
65
    if len(decBits) > 1:
66
        decDeg = float(decBits[0])
67
        if decBits[0].find("-") != -1:
68
            if len(decBits) > 1:
69
                decDeg = decDeg-(float(decBits[1])/60.0)
70
            if len(decBits) > 2:
71
                decDeg = decDeg-(float(decBits[2])/3600.0)
72
        else:
73
            if len(decBits) > 1:
74
                decDeg = decDeg+(float(decBits[1])/60.0)
75
            if len(decBits) > 2:
76
                decDeg = decDeg+(float(decBits[2])/3600.0)
77
    else:
78
        decDeg = float(decString)
79
80
    return decDeg
81
82
#-----------------------------------------------------------------------------
83
def decimal2hms(RADeg, delimiter):
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
    #if hours < 10 and hours >= 1:
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
    #if mins<10 and mins>=1:
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
    #if secs < 10 and secs>0.001:
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"        # So all to 3dp
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
        sDeg = str(int(sDeg)+1)
133
134
    return sHours+delimiter+sMins+delimiter+sSecs
135
136
#------------------------------------------------------------------------------
137
def decimal2dms(decDeg, delimiter):
138
    """Converts decimal degrees to string in Degrees:Minutes:Seconds format
139
    with user specified delimiter.
140

141
    @type decDeg: float
142
    @param decDeg: coordinate in decimal degrees
143
    @type delimiter: string
144
    @param delimiter: delimiter character in returned string
145
    @rtype: string
146
    @return: coordinate string in D:M:S format
147

148
    """
149
    # Positive
150
    if decDeg > 0:
151
        #if decDeg < 10 and decDeg>=1:
152
        if 1 <= decDeg < 10:
153
            sDeg = "0"+str(decDeg)[0]
154
        elif decDeg >= 10:
155
            sDeg = str(decDeg)[:2]
156
        elif decDeg < 1:
157
            sDeg = "00"
158
159
        if str(decDeg).find(".") == -1:
160
            mins = float(decDeg)*60.0
161
        else:
162
            mins = float(str(decDeg)[str(decDeg).index("."):])*60
163
        #if mins<10 and mins>=1:
164
        if 1 <= mins < 10:
165
            sMins = "0"+str(mins)[:1]
166
        elif mins >= 10:
167
            sMins = str(mins)[:2]
168
        elif mins < 1:
169
            sMins = "00"
170
171
        secs = (decDeg-(float(sDeg)+float(sMins)/60.0))*3600.0
172
        #if secs<10 and secs>0:
173
        if 0 < secs < 10:
174
            sSecs = "0"+str(secs)[:str(secs).find(".")+3]
175
        elif secs < 0.001:
176
            sSecs = "00.00"
177
        else:
178
            sSecs = str(secs)[:str(secs).find(".")+3]
179
        if len(sSecs) < 5:
180
            sSecs = sSecs+"0"        # So all to 2dp
181
182
        if float(sSecs) == 60.00:
183
            sSecs = "00.00"
184
            sMins = str(int(sMins)+1)
185
        if int(sMins) == 60:
186
            sMins = "00"
187
            sDeg = str(int(sDeg)+1)
188
189
        return "+"+sDeg+delimiter+sMins+delimiter+sSecs
190
191
    else:
192
        #if decDeg>-10 and decDeg<=-1:
193
        if -10 < decDeg <= -1:
194
            sDeg = "-0"+str(decDeg)[1]
195
        elif decDeg <= -10:
196
            sDeg = str(decDeg)[:3]
197
        elif decDeg > -1:
198
            sDeg = "-00"
199
200
        if str(decDeg).find(".") == -1:
201
            mins = float(decDeg)*-60.0
202
        else:
203
            mins = float(str(decDeg)[str(decDeg).index("."):])*60
204
        #if mins<10 and mins>=1:
205
        if 1 <= mins < 10:
206
            sMins = "0"+str(mins)[:1]
207
        elif mins >= 10:
208
            sMins = str(mins)[:2]
209
        elif mins < 1:
210
            sMins = "00"
211
212
        secs = (decDeg-(float(sDeg)-float(sMins)/60.0))*3600.0
213
        #if secs>-10 and secs<0:
214
        # so don't get minus sign
215
        if -10 < secs < 0:
216
            sSecs = "0"+str(secs)[1:str(secs).find(".")+3]
217
        elif secs > -0.001:
218
            sSecs = "00.00"
219
        else:
220
            sSecs = str(secs)[1:str(secs).find(".")+3]
221
        if len(sSecs) < 5:
222
            sSecs = sSecs+"0"        # So all to 2dp
223
224
        if float(sSecs) == 60.00:
225
            sSecs = "00.00"
226
            sMins = str(int(sMins)+1)
227
        if int(sMins) == 60:
228
            sMins = "00"
229
            sDeg = str(int(sDeg)-1)
230
231
        return sDeg+delimiter+sMins+delimiter+sSecs
232
233
#-----------------------------------------------------------------------------
234
def calcAngSepDeg(RA1, dec1, RA2, dec2):
235
    """Calculates the angular separation of two positions on the sky (specified
236
    in decimal degrees) in decimal degrees, assuming a tangent plane projection
237
    (so separation has to be <90 deg). Note that RADeg2, decDeg2 can be numpy
238
    arrays.
239

240
    @type RADeg1: float
241
    @param RADeg1: R.A. in decimal degrees for position 1
242
    @type decDeg1: float
243
    @param decDeg1: dec. in decimal degrees for position 1
244
    @type RADeg2: float or numpy array
245
    @param RADeg2: R.A. in decimal degrees for position 2
246
    @type decDeg2: float or numpy array
247
    @param decDeg2: dec. in decimal degrees for position 2
248
    @rtype: float or numpy array, depending upon type of RADeg2, decDeg2
249
    @return: angular separation in decimal degrees
250

251
    ***** Updates added by Alessandro Nastasi - 26/06/2014 *****
252

253
    Now RA and DEC can be entered indifferently as decimal degrees or as 
254
    hh:mm:ss.s or hh mm ss.s (provided that they are passed as STRING
255
    in the latter two cases). An internal routine recognizes the correct format 
256
    and converts them into decimal degrees.
257
 
258
    In addition, the tangent plane approximation restriction (i.e., dist < 90 deg) 
259
    has been removed and the complete formula is now implemented. 
260
    Pythagoras approximation is applied only for dist < 1 arcsec.
261
    """
262
    # ** ORIGINAL CODE **
263
    #
264
    #cRA = numpy.radians(RADeg1)
265
    #cDec = numpy.radians(decDeg1)
266
267
    #gRA = numpy.radians(RADeg2)
268
    #gDec = numpy.radians(decDeg2)
269
270
    #dRA = cRA-gRA
271
    #dDec = gDec-cDec
272
    #cosC = ((numpy.sin(gDec)*numpy.sin(cDec)) +
273
        #(numpy.cos(gDec)*numpy.cos(cDec) * numpy.cos(gRA-cRA)))
274
    #x = (numpy.cos(cDec)*numpy.sin(gRA-cRA))/cosC
275
    #y = (((numpy.cos(gDec)*numpy.sin(cDec)) - (numpy.sin(gDec) *
276
        #numpy.cos(cDec)*numpy.cos(gRA-cRA)))/cosC)
277
    #r = numpy.degrees(numpy.sqrt(x*x+y*y))
278
279
    RA1 = str(RA1).strip()
280
    dec1 = str(dec1).strip()
281
    RA2 = str(RA2).strip()
282
    dec2 = str(dec2).strip()
283
    
284
    inp_param = [str(RA1), dec1, RA2, dec2]
285
    newinp = []
286
    for x in inp_param: 
287
      newinp.append(string.replace(x, ":", " "))
288
      
289
    if str(newinp[0]).find(' ') >= 0:
290
      RADeg1 = hms2decimal(newinp[0], ' ')
291
      decDeg1 = dms2decimal(newinp[1], ' ')
292
    else:
293
      RADeg1 = float(newinp[0])
294
      decDeg1 = float(newinp[1])
295
296
    if str(newinp[2]).find(' ') >= 0:
297
      RADeg2 = hms2decimal(newinp[2], ' ')
298
      decDeg2 = dms2decimal(newinp[3], ' ')
299
    else:
300
      RADeg2 = float(newinp[2])
301
      decDeg2 = float(newinp[3])
302
        
303
    cRA = numpy.radians(RADeg1)
304
    cDec = numpy.radians(decDeg1)
305
306
    gRA = numpy.radians(RADeg2)
307
    gDec = numpy.radians(decDeg2)
308
      
309
    x=numpy.cos(cRA)*numpy.cos(cDec)*numpy.cos(gRA)*numpy.cos(gDec)
310
    y=numpy.sin(cRA)*numpy.cos(cDec)*numpy.sin(gRA)*numpy.cos(gDec)
311
    z=numpy.sin(cDec)*numpy.sin(gDec)
312
313
    rad=numpy.degrees(numpy.arccos(round(x+y+z,10)))
314
315
    dRA = cRA-gRA
316
    dDec = gDec-cDec
317
    #cosC = ((numpy.sin(gDec)*numpy.sin(cDec)) +
318
        #(numpy.cos(gDec)*numpy.cos(cDec) * numpy.cos(gRA-cRA)))
319
    #x = (numpy.cos(cDec)*numpy.sin(gRA-cRA))/cosC
320
    #y = (((numpy.cos(gDec)*numpy.sin(cDec)) - (numpy.sin(gDec) *
321
        #numpy.cos(cDec)*numpy.cos(gRA-cRA)))/cosC)
322
    #r = numpy.degrees(numpy.sqrt(x*x+y*y))
323
    
324
    # use Pythagoras approximation if rad < 1 arcsec
325
    if rad<0.000004848:
326
      rad=numpy.degrees(numpy.sqrt((numpy.cos(cDec)*dRA)**2+dDec**2))
327
328
    return rad
329
    
330
#-----------------------------------------------------------------------------
331
def calcAngSepDegPythagoras(RADeg1, decDeg1, RADeg2, decDeg2):
332
    # ** ORIGINAL CODE **
333
    #
334
    cRA = numpy.radians(RADeg1)
335
    cDec = numpy.radians(decDeg1)
336
337
    gRA = numpy.radians(RADeg2)
338
    gDec = numpy.radians(decDeg2)
339
340
    dRA = cRA-gRA
341
    dDec = gDec-cDec
342
    cosC = ((numpy.sin(gDec)*numpy.sin(cDec)) +
343
        (numpy.cos(gDec)*numpy.cos(cDec) * numpy.cos(gRA-cRA)))
344
    x = (numpy.cos(cDec)*numpy.sin(gRA-cRA))/cosC
345
    y = (((numpy.cos(gDec)*numpy.sin(cDec)) - (numpy.sin(gDec) *
346
        numpy.cos(cDec)*numpy.cos(gRA-cRA)))/cosC)
347
    r = numpy.degrees(numpy.sqrt(x*x+y*y))
348
  
349
    return r
350
351
#-----------------------------------------------------------------------------
352
def shiftRADec(ra1, dec1, deltaRA, deltaDec):
353
    """Computes new right ascension and declination shifted from the original
354
    by some delta RA and delta DEC. Input position is decimal degrees. Shifts
355
    (deltaRA, deltaDec) are arcseconds, and output is decimal degrees. Based on
356
    an IDL routine of the same name.
357

358
    @param ra1: float
359
    @type ra1: R.A. in decimal degrees
360
    @param dec1: float
361
    @type dec1: dec. in decimal degrees
362
    @param deltaRA: float
363
    @type deltaRA: shift in R.A. in arcseconds
364
    @param deltaDec: float
365
    @type deltaDec: shift in dec. in arcseconds
366
    @rtype: float [newRA, newDec]
367
    @return: shifted R.A. and dec.
368

369
    """
370
371
    d2r = math.pi/180.
372
    as2r = math.pi/648000.
373
374
    # Convert everything to radians
375
    rara1 = ra1*d2r
376
    dcrad1 = dec1*d2r
377
    shiftRArad = deltaRA*as2r
378
    shiftDCrad = deltaDec*as2r
379
380
    # Shift!
381
    deldec2 = 0.0
382
    sindis = math.sin(shiftRArad / 2.0)
383
    sindelRA = sindis / math.cos(dcrad1)
384
    delra = 2.0*math.asin(sindelRA) / d2r
385
386
    # Make changes
387
    ra2 = ra1+delra
388
    dec2 = dec1 +deltaDec / 3600.0
389
390
    return ra2, dec2
391
392
#-----------------------------------------------------------------------------
393
def convertCoords(inputSystem, outputSystem, coordX, coordY, epoch):
394
    """Converts specified coordinates (given in decimal degrees) between J2000,
395
    B1950, and Galactic.
396

397
    @type inputSystem: string
398
    @param inputSystem: system of the input coordinates (either "J2000",
399
        "B1950" or "GALACTIC")
400
    @type outputSystem: string
401
    @param outputSystem: system of the returned coordinates (either "J2000",
402
        "B1950" or "GALACTIC")
403
    @type coordX: float
404
    @param coordX: longitude coordinate in decimal degrees, e.g. R. A.
405
    @type coordY: float
406
    @param coordY: latitude coordinate in decimal degrees, e.g. dec.
407
    @type epoch: float
408
    @param epoch: epoch of the input coordinates
409
    @rtype: list
410
    @return: coordinates in decimal degrees in requested output system
411

412
    """
413
414
    if inputSystem=="J2000" or inputSystem=="B1950" or inputSystem=="GALACTIC":
415
        if outputSystem=="J2000" or outputSystem=="B1950" or \
416
            outputSystem=="GALACTIC":
417
418
            outCoords=wcscon.wcscon(wcscon.wcscsys(inputSystem),
419
                wcscon.wcscsys(outputSystem), 0, 0, coordX, coordY, epoch)
420
421
            return outCoords
422
423
    raise Exception("inputSystem and outputSystem must be 'J2000', 'B1950'"
424
                    "or 'GALACTIC'")
425
426
#-----------------------------------------------------------------------------
427
def calcRADecSearchBox(RADeg, decDeg, radiusSkyDeg):
428
    """Calculates minimum and maximum RA, dec coords needed to define a box
429
    enclosing a circle of radius radiusSkyDeg around the given RADeg, decDeg
430
    coordinates. Useful for freeform queries of e.g. SDSS, UKIDSS etc.. Uses
431
    L{calcAngSepDeg}, so has the same limitations.
432

433
    @type RADeg: float
434
    @param RADeg: RA coordinate of centre of search region
435
    @type decDeg: float
436
    @param decDeg: dec coordinate of centre of search region
437
    @type radiusSkyDeg: float
438
    @param radiusSkyDeg: radius in degrees on the sky used to define search
439
        region
440
    @rtype: list
441
    @return: [RAMin, RAMax, decMin, decMax] - coordinates in decimal degrees
442
        defining search box
443

444
    """
445
446
    tolerance = 1e-5  # in degrees on sky
447
    targetHalfSizeSkyDeg = radiusSkyDeg
448
    funcCalls = ["calcAngSepDeg(RADeg, decDeg, guess, decDeg)",
449
               "calcAngSepDeg(RADeg, decDeg, guess, decDeg)",
450
               "calcAngSepDeg(RADeg, decDeg, RADeg, guess)",
451
               "calcAngSepDeg(RADeg, decDeg, RADeg, guess)"]
452
    coords = [RADeg, RADeg, decDeg, decDeg]
453
    signs = [1.0, -1.0, 1.0, -1.0]
454
    results = []
455
    for f, c, sign in zip(funcCalls, coords, signs):
456
        # Initial guess range
457
        maxGuess = sign*targetHalfSizeSkyDeg*10.0
458
        minGuess = sign*targetHalfSizeSkyDeg/10.0
459
        guessStep = (maxGuess-minGuess)/10.0
460
        guesses = numpy.arange(minGuess+c, maxGuess+c, guessStep)
461
        for i in range(50):
462
            minSizeDiff = 1e6
463
            bestGuess = None
464
            for guess in guesses:
465
                sizeDiff = abs(eval(f)-targetHalfSizeSkyDeg)
466
                if sizeDiff < minSizeDiff:
467
                    minSizeDiff = sizeDiff
468
                    bestGuess = guess
469
            if minSizeDiff < tolerance:
470
                break
471
            else:
472
                guessRange = abs((maxGuess-minGuess))
473
                maxGuess = bestGuess+guessRange/4.0
474
                minGuess = bestGuess-guessRange/4.0
475
                guessStep = (maxGuess-minGuess)/10.0
476
                guesses = numpy.arange(minGuess, maxGuess, guessStep)
477
        results.append(bestGuess)
478
479
    RAMax = results[0]
480
    RAMin = results[1]
481
    decMax = results[2]
482
    decMin = results[3]
483
484
    return [RAMin, RAMax, decMin, decMax]