1"""Module for coordinate manipulation (conversions, calculations etc.). 2 3(c) 2007-2012 Matt Hilton 4 5(c) 2013-2016 Matt Hilton & Steven Boada 6 7""" 8 9import sys 10import math 11import numpy 12from PyWCSTools import wcscon 13#import IPython 14 15#----------------------------------------------------------------------------- 16def hms2decimal(RAString, delimiter): 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 # is it in HH:MM:SS format? 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#----------------------------------------------------------------------------- 46def dms2decimal(decString, delimiter): 47 """Converts a delimited string of Degrees:Minutes:Seconds format into 48 decimal degrees. 49 50 @type decString: string 51 @param decString: coordinate string in D:M:S format 52 @type delimiter: string 53 @param delimiter: delimiter character in decString 54 @rtype: float 55 @return: coordinate in decimal degrees 56 57 """ 58 # is it in DD:MM:SS format? 59 if delimiter == "": 60 decBits = str(decString).split() 61 else: 62 decBits = str(decString).split(delimiter) 63 if len(decBits) > 1: 64 decDeg = float(decBits[0]) 65 if decBits[0].find("-") != -1: 66 if len(decBits) > 1: 67 decDeg = decDeg-(float(decBits[1])/60.0) 68 if len(decBits) > 2: 69 decDeg = decDeg-(float(decBits[2])/3600.0) 70 else: 71 if len(decBits) > 1: 72 decDeg = decDeg+(float(decBits[1])/60.0) 73 if len(decBits) > 2: 74 decDeg = decDeg+(float(decBits[2])/3600.0) 75 else: 76 decDeg = float(decString) 77 78 return decDeg 79 80#----------------------------------------------------------------------------- 81def decimal2hms(RADeg, delimiter): 82 """Converts decimal degrees to string in Hours:Minutes:Seconds format with 83 user specified delimiter. 84 85 @type RADeg: float 86 @param RADeg: coordinate in decimal degrees 87 @type delimiter: string 88 @param delimiter: delimiter character in returned string 89 @rtype: string 90 @return: coordinate string in H:M:S format 91 92 """ 93 hours = (RADeg/360.0)*24 94 #if hours < 10 and hours >= 1: 95 if 1 <= hours < 10: 96 sHours = "0"+str(hours)[0] 97 elif hours >= 10: 98 sHours = str(hours)[:2] 99 elif hours < 1: 100 sHours = "00" 101 102 if str(hours).find(".") == -1: 103 mins = float(hours)*60.0 104 else: 105 mins = float(str(hours)[str(hours).index("."):])*60.0 106 #if mins<10 and mins>=1: 107 if 1 <= mins<10: 108 sMins = "0"+str(mins)[:1] 109 elif mins >= 10: 110 sMins = str(mins)[:2] 111 elif mins < 1: 112 sMins = "00" 113 114 secs = (hours-(float(sHours)+float(sMins)/60.0))*3600.0 115 #if secs < 10 and secs>0.001: 116 if 0.001 < secs < 10: 117 sSecs = "0"+str(secs)[:str(secs).find(".")+4] 118 elif secs < 0.0001: 119 sSecs = "00.000" 120 else: 121 sSecs = str(secs)[:str(secs).find(".")+4] 122 if len(sSecs) < 5: 123 sSecs = sSecs+"00" # So all to 3dp 124 125 if float(sSecs) == 60.000: 126 sSecs = "00.00" 127 sMins = str(int(sMins)+1) 128 if int(sMins) == 60: 129 sMins = "00" 130 sHours = str(int(sHours)+1) 131 132 return sHours+delimiter+sMins+delimiter+sSecs 133 134#------------------------------------------------------------------------------ 135def decimal2dms(decDeg, delimiter): 136 """Converts decimal degrees to string in Degrees:Minutes:Seconds format 137 with user specified delimiter. 138 139 @type decDeg: float 140 @param decDeg: coordinate in decimal degrees 141 @type delimiter: string 142 @param delimiter: delimiter character in returned string 143 @rtype: string 144 @return: coordinate string in D:M:S format 145 146 """ 147 # Positive 148 if decDeg > 0: 149 #if decDeg < 10 and decDeg>=1: 150 if 1 <= decDeg < 10: 151 sDeg = "0"+str(decDeg)[0] 152 elif decDeg >= 10: 153 sDeg = str(decDeg)[:2] 154 elif decDeg < 1: 155 sDeg = "00" 156 157 if str(decDeg).find(".") == -1: 158 mins = float(decDeg)*60.0 159 else: 160 mins = float(str(decDeg)[str(decDeg).index("."):])*60 161 #if mins<10 and mins>=1: 162 if 1 <= mins < 10: 163 sMins = "0"+str(mins)[:1] 164 elif mins >= 10: 165 sMins = str(mins)[:2] 166 elif mins < 1: 167 sMins = "00" 168 169 secs = (decDeg-(float(sDeg)+float(sMins)/60.0))*3600.0 170 #if secs<10 and secs>0: 171 if 0 < secs < 10: 172 sSecs = "0"+str(secs)[:str(secs).find(".")+3] 173 elif secs < 0.001: 174 sSecs = "00.00" 175 else: 176 sSecs = str(secs)[:str(secs).find(".")+3] 177 if len(sSecs) < 5: 178 sSecs = sSecs+"0" # So all to 2dp 179 180 if float(sSecs) == 60.00: 181 sSecs = "00.00" 182 sMins = str(int(sMins)+1) 183 if int(sMins) == 60: 184 sMins = "00" 185 sDeg = str(int(sDeg)+1) 186 187 return "+"+sDeg+delimiter+sMins+delimiter+sSecs 188 189 else: 190 #if decDeg>-10 and decDeg<=-1: 191 if -10 < decDeg <= -1: 192 sDeg = "-0"+str(decDeg)[1] 193 elif decDeg <= -10: 194 sDeg = str(decDeg)[:3] 195 elif decDeg > -1: 196 sDeg = "-00" 197 198 if str(decDeg).find(".") == -1: 199 mins = float(decDeg)*-60.0 200 else: 201 mins = float(str(decDeg)[str(decDeg).index("."):])*60 202 #if mins<10 and mins>=1: 203 if 1 <= mins < 10: 204 sMins = "0"+str(mins)[:1] 205 elif mins >= 10: 206 sMins = str(mins)[:2] 207 elif mins < 1: 208 sMins = "00" 209 210 secs = (decDeg-(float(sDeg)-float(sMins)/60.0))*3600.0 211 #if secs>-10 and secs<0: 212 # so don't get minus sign 213 if -10 < secs < 0: 214 sSecs = "0"+str(secs)[1:str(secs).find(".")+3] 215 elif secs > -0.001: 216 sSecs = "00.00" 217 else: 218 sSecs = str(secs)[1:str(secs).find(".")+3] 219 if len(sSecs) < 5: 220 sSecs = sSecs+"0" # So all to 2dp 221 222 if float(sSecs) == 60.00: 223 sSecs = "00.00" 224 sMins = str(int(sMins)+1) 225 if int(sMins) == 60: 226 sMins = "00" 227 sDeg = str(int(sDeg)-1) 228 229 return sDeg+delimiter+sMins+delimiter+sSecs 230 231#----------------------------------------------------------------------------- 232def calcAngSepDeg(RADeg1, decDeg1, RADeg2, decDeg2): 233 """Calculates the angular separation of two positions on the sky (specified 234 in decimal degrees) in decimal degrees. Note that RADeg2, decDeg2 can be numpy 235 arrays. 236 237 @type RADeg1: float 238 @param RADeg1: R.A. in decimal degrees for position 1 239 @type decDeg1: float 240 @param decDeg1: dec. in decimal degrees for position 1 241 @type RADeg2: float or numpy array 242 @param RADeg2: R.A. in decimal degrees for position 2 243 @type decDeg2: float or numpy array 244 @param decDeg2: dec. in decimal degrees for position 2 245 @rtype: float or numpy array, depending upon type of RADeg2, decDeg2 246 @return: angular separation in decimal degrees 247 248 """ 249 250 a=numpy.sin(numpy.radians(decDeg1))*numpy.sin(numpy.radians(decDeg2))+numpy.cos(numpy.radians(decDeg1))*numpy.cos(numpy.radians(decDeg2))*numpy.cos(numpy.radians(RADeg1-RADeg2)) 251 mask=numpy.greater(a, 1.0) 252 if mask.sum() > 0: 253 if type(a) == numpy.ndarray: 254 a[mask]=1.0 255 else: 256 a=1.0 257 mask=numpy.less(a, -1.0) 258 if mask.sum() > 0: 259 if type(a) == numpy.ndarray: 260 a[mask]=-1.0 261 else: 262 a=-1.0 263 r=numpy.degrees(numpy.arccos(a)) 264 265 # Above gives nan when RADeg1, decDeg1 == RADeg1, decDeg2 266 indexList=numpy.where(numpy.isnan(r) == True)[0] 267 tolerance=1e-6 268 if len(indexList) > 0: 269 for index in indexList: 270 if type(r) == numpy.ndarray: 271 if type(RADeg2) == numpy.ndarray: 272 if abs(RADeg1 - RADeg2[index]) < tolerance and abs(decDeg1 -decDeg2[index]) < tolerance: 273 r[index]=0.0 274 else: 275 raise Exception("astCoords: calcAngSepDeg - encountered nan not due to equal RADeg, decDeg coords") 276 elif type(RADeg1) == numpy.ndarray: 277 if abs(RADeg2 - RADeg1[index]) < tolerance and abs(decDeg2 -decDeg1[index]) < tolerance: 278 r[index]=0.0 279 else: 280 raise Exception("astCoords: calcAngSepDeg - encountered nan not due to equal RADeg, decDeg coords") 281 else: 282 r=0.0 283 284 return r 285 286#----------------------------------------------------------------------------- 287def shiftRADec(ra1, dec1, deltaRA, deltaDec): 288 """Computes new right ascension and declination shifted from the original 289 by some delta RA and delta DEC. Input position is decimal degrees. Shifts 290 (deltaRA, deltaDec) are arcseconds, and output is decimal degrees. Based on 291 an IDL routine of the same name. 292 293 @param ra1: float 294 @type ra1: R.A. in decimal degrees 295 @param dec1: float 296 @type dec1: dec. in decimal degrees 297 @param deltaRA: float 298 @type deltaRA: shift in R.A. in arcseconds 299 @param deltaDec: float 300 @type deltaDec: shift in dec. in arcseconds 301 @rtype: float [newRA, newDec] 302 @return: shifted R.A. and dec. 303 304 """ 305 306 d2r = math.pi/180. 307 as2r = math.pi/648000. 308 309 # Convert everything to radians 310 rara1 = ra1*d2r 311 dcrad1 = dec1*d2r 312 shiftRArad = deltaRA*as2r 313 shiftDCrad = deltaDec*as2r 314 315 # Shift! 316 deldec2 = 0.0 317 sindis = math.sin(shiftRArad / 2.0) 318 sindelRA = sindis / math.cos(dcrad1) 319 delra = 2.0*math.asin(sindelRA) / d2r 320 321 # Make changes 322 ra2 = ra1+delra 323 dec2 = dec1 +deltaDec / 3600.0 324 325 return ra2, dec2 326 327#----------------------------------------------------------------------------- 328def convertCoords(inputSystem, outputSystem, coordX, coordY, epoch): 329 """Converts specified coordinates (given in decimal degrees) between J2000, 330 B1950, and Galactic. 331 332 @type inputSystem: string 333 @param inputSystem: system of the input coordinates (either "J2000", 334 "B1950" or "GALACTIC") 335 @type outputSystem: string 336 @param outputSystem: system of the returned coordinates (either "J2000", 337 "B1950" or "GALACTIC") 338 @type coordX: float 339 @param coordX: longitude coordinate in decimal degrees, e.g. R. A. 340 @type coordY: float 341 @param coordY: latitude coordinate in decimal degrees, e.g. dec. 342 @type epoch: float 343 @param epoch: epoch of the input coordinates 344 @rtype: list 345 @return: coordinates in decimal degrees in requested output system 346 347 """ 348 349 if inputSystem=="J2000" or inputSystem=="B1950" or inputSystem=="GALACTIC": 350 if outputSystem=="J2000" or outputSystem=="B1950" or \ 351 outputSystem=="GALACTIC": 352 353 outCoords=wcscon.wcscon(wcscon.wcscsys(inputSystem), 354 wcscon.wcscsys(outputSystem), 0, 0, float(coordX), float(coordY), epoch) 355 356 return outCoords 357 358 raise Exception("inputSystem and outputSystem must be 'J2000', 'B1950'" 359 "or 'GALACTIC'") 360 361#----------------------------------------------------------------------------- 362def calcRADecSearchBox(RADeg, decDeg, radiusSkyDeg): 363 """Calculates minimum and maximum RA, dec coords needed to define a box 364 enclosing a circle of radius radiusSkyDeg around the given RADeg, decDeg 365 coordinates. Useful for freeform queries of e.g. SDSS, UKIDSS etc.. Uses 366 L{calcAngSepDeg}, so has the same limitations. 367 368 @type RADeg: float 369 @param RADeg: RA coordinate of centre of search region 370 @type decDeg: float 371 @param decDeg: dec coordinate of centre of search region 372 @type radiusSkyDeg: float 373 @param radiusSkyDeg: radius in degrees on the sky used to define search 374 region 375 @rtype: list 376 @return: [RAMin, RAMax, decMin, decMax] - coordinates in decimal degrees 377 defining search box 378 379 """ 380 381 tolerance = 1e-5 # in degrees on sky 382 targetHalfSizeSkyDeg = radiusSkyDeg 383 funcCalls = ["calcAngSepDeg(RADeg, decDeg, guess, decDeg)", 384 "calcAngSepDeg(RADeg, decDeg, guess, decDeg)", 385 "calcAngSepDeg(RADeg, decDeg, RADeg, guess)", 386 "calcAngSepDeg(RADeg, decDeg, RADeg, guess)"] 387 coords = [RADeg, RADeg, decDeg, decDeg] 388 signs = [1.0, -1.0, 1.0, -1.0] 389 results = [] 390 for f, c, sign in zip(funcCalls, coords, signs): 391 # Initial guess range 392 maxGuess = sign*targetHalfSizeSkyDeg*80.0 393 minGuess = sign*targetHalfSizeSkyDeg/80.0 394 #guessStep = (maxGuess-minGuess)/10.0 395 guesses = numpy.linspace(minGuess+c, maxGuess+c, 1000) 396 converged=False 397 for i in range(50): 398 minSizeDiff = 1e6 399 bestGuess = None 400 for guess in guesses: 401 sizeDiff = abs(eval(f)-targetHalfSizeSkyDeg) 402 if sizeDiff < minSizeDiff: 403 minSizeDiff = sizeDiff 404 bestGuess = guess 405 if minSizeDiff < tolerance: 406 converged=True 407 break 408 else: 409 #print sizeDiff, bestGuess, bestGuess-minGuess, bestGuess-maxGuess 410 if bestGuess == None: 411 raise Exception("bestGuess is None") 412 guessRange = abs((maxGuess-minGuess)) 413 maxGuess = bestGuess+guessRange/4.0 414 minGuess = bestGuess-guessRange/4.0 415 # Stop us from searching the wrong side of the coordinate 416 if sign == 1: 417 if minGuess < c: 418 minGuess=c+tolerance 419 if sign == -1: 420 if maxGuess > c: 421 maxGuess=c-tolerance 422 #guessStep = (maxGuess-minGuess)/20.0 423 guesses = numpy.linspace(minGuess, maxGuess, 1000) 424 if converged == False: 425 raise Exception("calcRADecSearchBox failed to converge") 426 results.append(bestGuess) 427 428 RAMax = results[0] 429 RAMin = results[1] 430 decMax = results[2] 431 decMin = results[3] 432 433 # Sanity check 434 if (RAMax-RAMin)+(2*tolerance) < 2*targetHalfSizeSkyDeg: 435 raise Exception("calcRADecSearchBox failed sanity check") 436 437 return [RAMin, RAMax, decMin, decMax] 438 439