Package astLib :: Module astCoords
[hide private]
[frames] | no frames]

Source Code for Module astLib.astCoords

  1  """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  from PyWCSTools import wcscon 
 15   
 16  #----------------------------------------------------------------------------- 
17 -def hms2decimal(RAString, delimiter):
18 """Converts a delimited string of Hours:Minutes:Seconds format into decimal 19 degrees. 20 21 @type RAString: string 22 @param RAString: coordinate string in H:M:S format 23 @type delimiter: string 24 @param delimiter: delimiter character in RAString 25 @rtype: float 26 @return: coordinate in decimal degrees 27 28 """ 29 # is it in HH:MM:SS format? 30 if delimiter == "": 31 RABits = str(RAString).split() 32 else: 33 RABits = str(RAString).split(delimiter) 34 if len(RABits) > 1: 35 RAHDecimal = float(RABits[0]) 36 if len(RABits) > 1: 37 RAHDecimal = RAHDecimal+(float(RABits[1])/60.0) 38 if len(RABits) > 2: 39 RAHDecimal = RAHDecimal+(float(RABits[2])/3600.0) 40 RADeg = (RAHDecimal/24.0)*360.0 41 else: 42 RADeg = float(RAString) 43 44 return RADeg
45 46 #-----------------------------------------------------------------------------
47 -def dms2decimal(decString, delimiter):
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 # is it in DD:MM:SS format? 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 -def decimal2hms(RADeg, delimiter):
83 """Converts decimal degrees to string in Hours:Minutes:Seconds format with 84 user specified delimiter. 85 86 @type RADeg: float 87 @param RADeg: coordinate in decimal degrees 88 @type delimiter: string 89 @param delimiter: delimiter character in returned string 90 @rtype: string 91 @return: coordinate string in H:M:S format 92 93 """ 94 hours = (RADeg/360.0)*24 95 #if hours < 10 and hours >= 1: 96 if 1 <= hours < 10: 97 sHours = "0"+str(hours)[0] 98 elif hours >= 10: 99 sHours = str(hours)[:2] 100 elif hours < 1: 101 sHours = "00" 102 103 if str(hours).find(".") == -1: 104 mins = float(hours)*60.0 105 else: 106 mins = float(str(hours)[str(hours).index("."):])*60.0 107 #if mins<10 and mins>=1: 108 if 1 <= mins<10: 109 sMins = "0"+str(mins)[:1] 110 elif mins >= 10: 111 sMins = str(mins)[:2] 112 elif mins < 1: 113 sMins = "00" 114 115 secs = (hours-(float(sHours)+float(sMins)/60.0))*3600.0 116 #if secs < 10 and secs>0.001: 117 if 0.001 < secs < 10: 118 sSecs = "0"+str(secs)[:str(secs).find(".")+4] 119 elif secs < 0.0001: 120 sSecs = "00.001" 121 else: 122 sSecs = str(secs)[:str(secs).find(".")+4] 123 if len(sSecs) < 5: 124 sSecs = sSecs+"00" # So all to 3dp 125 126 if float(sSecs) == 60.000: 127 sSecs = "00.00" 128 sMins = str(int(sMins)+1) 129 if int(sMins) == 60: 130 sMins = "00" 131 sDeg = str(int(sDeg)+1) 132 133 return sHours+delimiter+sMins+delimiter+sSecs
134 135 #------------------------------------------------------------------------------
136 -def decimal2dms(decDeg, delimiter):
137 """Converts decimal degrees to string in Degrees:Minutes:Seconds format 138 with user specified delimiter. 139 140 @type decDeg: float 141 @param decDeg: coordinate in decimal degrees 142 @type delimiter: string 143 @param delimiter: delimiter character in returned string 144 @rtype: string 145 @return: coordinate string in D:M:S format 146 147 """ 148 # Positive 149 if decDeg > 0: 150 #if decDeg < 10 and decDeg>=1: 151 if 1 <= decDeg < 10: 152 sDeg = "0"+str(decDeg)[0] 153 elif decDeg >= 10: 154 sDeg = str(decDeg)[:2] 155 elif decDeg < 1: 156 sDeg = "00" 157 158 if str(decDeg).find(".") == -1: 159 mins = float(decDeg)*60.0 160 else: 161 mins = float(str(decDeg)[str(decDeg).index("."):])*60 162 #if mins<10 and mins>=1: 163 if 1 <= mins < 10: 164 sMins = "0"+str(mins)[:1] 165 elif mins >= 10: 166 sMins = str(mins)[:2] 167 elif mins < 1: 168 sMins = "00" 169 170 secs = (decDeg-(float(sDeg)+float(sMins)/60.0))*3600.0 171 #if secs<10 and secs>0: 172 if 0 < secs < 10: 173 sSecs = "0"+str(secs)[:str(secs).find(".")+3] 174 elif secs < 0.001: 175 sSecs = "00.00" 176 else: 177 sSecs = str(secs)[:str(secs).find(".")+3] 178 if len(sSecs) < 5: 179 sSecs = sSecs+"0" # So all to 2dp 180 181 if float(sSecs) == 60.00: 182 sSecs = "00.00" 183 sMins = str(int(sMins)+1) 184 if int(sMins) == 60: 185 sMins = "00" 186 sDeg = str(int(sDeg)+1) 187 188 return "+"+sDeg+delimiter+sMins+delimiter+sSecs 189 190 else: 191 #if decDeg>-10 and decDeg<=-1: 192 if -10 < decDeg <= -1: 193 sDeg = "-0"+str(decDeg)[1] 194 elif decDeg <= -10: 195 sDeg = str(decDeg)[:3] 196 elif decDeg > -1: 197 sDeg = "-00" 198 199 if str(decDeg).find(".") == -1: 200 mins = float(decDeg)*-60.0 201 else: 202 mins = float(str(decDeg)[str(decDeg).index("."):])*60 203 #if mins<10 and mins>=1: 204 if 1 <= mins < 10: 205 sMins = "0"+str(mins)[:1] 206 elif mins >= 10: 207 sMins = str(mins)[:2] 208 elif mins < 1: 209 sMins = "00" 210 211 secs = (decDeg-(float(sDeg)-float(sMins)/60.0))*3600.0 212 #if secs>-10 and secs<0: 213 # so don't get minus sign 214 if -10 < secs < 0: 215 sSecs = "0"+str(secs)[1:str(secs).find(".")+3] 216 elif secs > -0.001: 217 sSecs = "00.00" 218 else: 219 sSecs = str(secs)[1:str(secs).find(".")+3] 220 if len(sSecs) < 5: 221 sSecs = sSecs+"0" # So all to 2dp 222 223 if float(sSecs) == 60.00: 224 sSecs = "00.00" 225 sMins = str(int(sMins)+1) 226 if int(sMins) == 60: 227 sMins = "00" 228 sDeg = str(int(sDeg)-1) 229 230 return sDeg+delimiter+sMins+delimiter+sSecs
231 232 #-----------------------------------------------------------------------------
233 -def calcAngSepDeg(RADeg1, decDeg1, RADeg2, decDeg2):
234 """Calculates the angular separation of two positions on the sky (specified 235 in decimal degrees) in decimal degrees, assuming a tangent plane projection 236 (so separation has to be <90 deg). Note that RADeg2, decDeg2 can be numpy 237 arrays. 238 239 @type RADeg1: float 240 @param RADeg1: R.A. in decimal degrees for position 1 241 @type decDeg1: float 242 @param decDeg1: dec. in decimal degrees for position 1 243 @type RADeg2: float or numpy array 244 @param RADeg2: R.A. in decimal degrees for position 2 245 @type decDeg2: float or numpy array 246 @param decDeg2: dec. in decimal degrees for position 2 247 @rtype: float or numpy array, depending upon type of RADeg2, decDeg2 248 @return: angular separation in decimal degrees 249 250 """ 251 cRA = numpy.radians(RADeg1) 252 cDec = numpy.radians(decDeg1) 253 254 gRA = numpy.radians(RADeg2) 255 gDec = numpy.radians(decDeg2) 256 257 dRA = cRA-gRA 258 dDec = gDec-cDec 259 cosC = ((numpy.sin(gDec)*numpy.sin(cDec)) + 260 (numpy.cos(gDec)*numpy.cos(cDec) * numpy.cos(gRA-cRA))) 261 x = (numpy.cos(cDec)*numpy.sin(gRA-cRA))/cosC 262 y = (((numpy.cos(gDec)*numpy.sin(cDec)) - (numpy.sin(gDec) * 263 numpy.cos(cDec)*numpy.cos(gRA-cRA)))/cosC) 264 r = numpy.degrees(numpy.sqrt(x*x+y*y)) 265 266 return r
267 268 #-----------------------------------------------------------------------------
269 -def shiftRADec(ra1, dec1, deltaRA, deltaDec):
270 """Computes new right ascension and declination shifted from the original 271 by some delta RA and delta DEC. Input position is decimal degrees. Shifts 272 (deltaRA, deltaDec) are arcseconds, and output is decimal degrees. Based on 273 an IDL routine of the same name. 274 275 @param ra1: float 276 @type ra1: R.A. in decimal degrees 277 @param dec1: float 278 @type dec1: dec. in decimal degrees 279 @param deltaRA: float 280 @type deltaRA: shift in R.A. in arcseconds 281 @param deltaDec: float 282 @type deltaDec: shift in dec. in arcseconds 283 @rtype: float [newRA, newDec] 284 @return: shifted R.A. and dec. 285 286 """ 287 288 d2r = math.pi/180. 289 as2r = math.pi/648000. 290 291 # Convert everything to radians 292 rara1 = ra1*d2r 293 dcrad1 = dec1*d2r 294 shiftRArad = deltaRA*as2r 295 shiftDCrad = deltaDec*as2r 296 297 # Shift! 298 deldec2 = 0.0 299 sindis = math.sin(shiftRArad / 2.0) 300 sindelRA = sindis / math.cos(dcrad1) 301 delra = 2.0*math.asin(sindelRA) / d2r 302 303 # Make changes 304 ra2 = ra1+delra 305 dec2 = dec1 +deltaDec / 3600.0 306 307 return ra2, dec2
308 309 #-----------------------------------------------------------------------------
310 -def convertCoords(inputSystem, outputSystem, coordX, coordY, epoch):
311 """Converts specified coordinates (given in decimal degrees) between J2000, 312 B1950, and Galactic. 313 314 @type inputSystem: string 315 @param inputSystem: system of the input coordinates (either "J2000", 316 "B1950" or "GALACTIC") 317 @type outputSystem: string 318 @param outputSystem: system of the returned coordinates (either "J2000", 319 "B1950" or "GALACTIC") 320 @type coordX: float 321 @param coordX: longitude coordinate in decimal degrees, e.g. R. A. 322 @type coordY: float 323 @param coordY: latitude coordinate in decimal degrees, e.g. dec. 324 @type epoch: float 325 @param epoch: epoch of the input coordinates 326 @rtype: list 327 @return: coordinates in decimal degrees in requested output system 328 329 """ 330 331 if inputSystem=="J2000" or inputSystem=="B1950" or inputSystem=="GALACTIC": 332 if outputSystem=="J2000" or outputSystem=="B1950" or \ 333 outputSystem=="GALACTIC": 334 335 outCoords=wcscon.wcscon(wcscon.wcscsys(inputSystem), 336 wcscon.wcscsys(outputSystem), 0, 0, coordX, coordY, epoch) 337 338 return outCoords 339 340 raise Exception("inputSystem and outputSystem must be 'J2000', 'B1950'" 341 "or 'GALACTIC'")
342 343 #-----------------------------------------------------------------------------
344 -def calcRADecSearchBox(RADeg, decDeg, radiusSkyDeg):
345 """Calculates minimum and maximum RA, dec coords needed to define a box 346 enclosing a circle of radius radiusSkyDeg around the given RADeg, decDeg 347 coordinates. Useful for freeform queries of e.g. SDSS, UKIDSS etc.. Uses 348 L{calcAngSepDeg}, so has the same limitations. 349 350 @type RADeg: float 351 @param RADeg: RA coordinate of centre of search region 352 @type decDeg: float 353 @param decDeg: dec coordinate of centre of search region 354 @type radiusSkyDeg: float 355 @param radiusSkyDeg: radius in degrees on the sky used to define search 356 region 357 @rtype: list 358 @return: [RAMin, RAMax, decMin, decMax] - coordinates in decimal degrees 359 defining search box 360 361 """ 362 363 tolerance = 1e-5 # in degrees on sky 364 targetHalfSizeSkyDeg = radiusSkyDeg 365 funcCalls = ["calcAngSepDeg(RADeg, decDeg, guess, decDeg)", 366 "calcAngSepDeg(RADeg, decDeg, guess, decDeg)", 367 "calcAngSepDeg(RADeg, decDeg, RADeg, guess)", 368 "calcAngSepDeg(RADeg, decDeg, RADeg, guess)"] 369 coords = [RADeg, RADeg, decDeg, decDeg] 370 signs = [1.0, -1.0, 1.0, -1.0] 371 results = [] 372 for f, c, sign in zip(funcCalls, coords, signs): 373 # Initial guess range 374 maxGuess = sign*targetHalfSizeSkyDeg*10.0 375 minGuess = sign*targetHalfSizeSkyDeg/10.0 376 guessStep = (maxGuess-minGuess)/10.0 377 guesses = numpy.arange(minGuess+c, maxGuess+c, guessStep) 378 for i in range(50): 379 minSizeDiff = 1e6 380 bestGuess = None 381 for guess in guesses: 382 sizeDiff = abs(eval(f)-targetHalfSizeSkyDeg) 383 if sizeDiff < minSizeDiff: 384 minSizeDiff = sizeDiff 385 bestGuess = guess 386 if minSizeDiff < tolerance: 387 break 388 else: 389 guessRange = abs((maxGuess-minGuess)) 390 maxGuess = bestGuess+guessRange/4.0 391 minGuess = bestGuess-guessRange/4.0 392 guessStep = (maxGuess-minGuess)/10.0 393 guesses = numpy.arange(minGuess, maxGuess, guessStep) 394 results.append(bestGuess) 395 396 RAMax = results[0] 397 RAMin = results[1] 398 decMax = results[2] 399 decMin = results[3] 400 401 return [RAMin, RAMax, decMin, decMax]
402