Right ascension and Declination#

Right ascension is often given in hours-minutes-seconds (HMS) notation, because it was convenient to calculate when a star would appear over the horizon. A full circle in HMS notation is 24 hours, which means 1 hour in HMS notation is equal to 15 degrees.

Each hour is split into 60 minutes and each minute into 60 seconds.

You can convert 23 hours, 12 minutes and 6 seconds (written as 23:12:06 or 23h12m06s) to degrees like this:

>>> print(15*(23 + 12/60 + 6/(60*60)))
348.025

print(15*(23 + 12/60 + 6/(60*60)))

348.025

Declination, on the other hand, is traditionally recorded in degrees-minutes-seconds (DMS) notation. A full circle is 360 degrees, each degree has 60 arcminutes and each arcminute has 60 arcseconds.

For example: 73 degrees, 21 arcminutes and 14.4 arcseconds (written 73:21:14.4 or 73° 21’ 14.4” or 73d21m14.4s) can be converted to decimal degrees like this:

>>> print(73 + 21/60 + 14.4/(60*60))
73.354

With negative angles like -5° 31’ 12” the negation applies to the whole angle, including arcminutes and arcseconds:

>>> print(-1*(5 + 31/60 + 12/(60*60)))
-5.52

Note: arcminutes and minutes are different! The arcminutes and arcseconds in DMS are not the same as the minutes and seconds in HMS! A minute in HMS is equal to 15 arcminutes in DMS and a second is equal to 15 arcseconds.

Write two functions, one that converts right ascension from HMS to decimal degrees, called hms2dec, and another that converts declination from DMS to decimal degrees, called dms2dec.

Right ascension is always an angle from 0 to 24 hours and declination is always an angle from -90° to +90°.

Your hms2dec function should work like this:

>>> hms2dec(23, 12, 6)
348.025

And your dms2dec function should work like this:

>>> dms2dec(22, 57, 18)
22.955

It should also work with negative angles:

>>> dms2dec(-66, 5, 5.1)
-66.08475

# Write your hms2dec and dms2dec functions here
def hms2dec(h: float, m:float, s:float) -> float:
    pass

def dms2dec(d: float, m: float, s:float) -> float:
    pass

# You can use this to test your function.
# The first example from the question
print(hms2dec(23, 12, 6))

# The second example from the question
print(dms2dec(22, 57, 18))

# The third example from the question
print(dms2dec(-66, 5, 5.1))

348.025
22.955
-66.08475

Before we can crossmatch our two catalogues we first have to import the raw data. Every astronomy catalogue tends to have its own unique format so we’ll need to look at how to do this with each one individually.

We’ll look at the AT20G bright source sample survey first. The raw data we’ll be using is the file table2.dat from this page in the VizieR archives, but we’ll use the filename bss.dat from now on.

Every catalogue in VizieR has a detailed README file that gives you the exact format of each table in the catalogue.

The full catalogue of bright radio sources contains 320 objects. The first few rows look like this (scroll right to see it all):

 1  00 04 35.65 -47 36 19.1   0.87 0.04 0.97 0.06  0.90 0.04                0.995 0.030            17.63 Q 1.F.11.C  PKS 0002-478
 2  00 10 35.92 -30 27 48.3   0.74 0.03 0.72 0.04  0.63 0.03  0.315 0.009   0.419 0.013 1.19  La01 19.59 Q 1.F.11..  PKS 0008-307
 3* 00 11 01.27 -26 12 33.1   0.64 0.07 0.82 0.07  0.69 0.03  0.210 0.006               1.096 Wr83 19.53 Q 4.F.44.C  PKS 0008-264

The catalogue is organised in fixed-width columns, with the format of the columns being:

1: Object catalogue ID number (sometimes with an asterisk)
2-4: Right ascension in HMS notation
5-7: Declination in DMS notation
8-: Other information, including spectral intensities

We only need coordinates for crossmatching. We can load specific columns with the usecols argument in NumPy’s loadtxt function:

import numpy as np
cat = np.loadtxt('bss.dat', usecols=range(1, 7))
print(cat[0])

We’ve skipped the ID column, since the ID number is always the same as the row number.

Fixed-width columns and loadtxt

loadtxt does not work for fixed-width columns if values are missing. Since there are no missing ID, RA and dec values it is fine for loading the first few columns of the BSS catalogue.

# Write your crossmatch function here.
from scipy.spatial import cKDTree
import numpy as np

def sign(d,m,s):
  ret = -1 if d<0 or m<0 or s<0 else +1
  return ret

def hms2dec(h,m,s):
  return 15*(h + m/60. + s/3600.)

def dms2dec(d,m,s):
  return sign(d,m,s)*(abs(d) + abs(m)/60. + abs(s)/3600.)

def import_bss():
  dat = np.loadtxt('../data/bss.dat', usecols=range(1,7))
  ret = []
  for idx, row in enumerate(dat):
    hRA, mRA, sRA, hDEC, mDEC, sDEC = row
    RA = hms2dec(hRA, mRA, sRA)
    DEC = dms2dec(hDEC, mDEC, sDEC)
    ret.append( (idx+1, RA, DEC) )

  return ret

def import_super():
  dat = np.loadtxt('../data/super.csv', delimiter=',', skiprows=1, usecols=[0,1])
  ret = zip(np.linspace(1,len(dat), len(dat)), dat[:,0], dat[:,1])
  return list(ret)

def angular_dist(ra1, dec1, ra2, dec2):
  dalpha = np.abs(np.radians(ra1-ra2))
  delta1 = np.radians(dec1)
  delta2 = np.radians(dec2)

  a = np.sin(0.5*np.abs(delta1-delta2))**2

  b = np.cos(delta1)*np.cos(delta2)*(np.sin(0.5*dalpha)**2)

  d = 2*np.arcsin(np.sqrt(a+b))

  return np.degrees(d)

def crossmatch(cat1, cat2, max_dist):
  dat1 = np.array(cat1)[:,1:]
  dat2 = np.array(cat2)[:,1:]
  tree = cKDTree(dat2)
  d2d, idx = tree.query(dat1, distance_upper_bound=max_dist*4) # allow 400% margin of error

  dists = []
  for id1, id2 in enumerate(idx):
    if id2<len(dat2):
      ra1, dec1 = dat1[id1]
      ra2, dec2 = dat2[id2]
      dist = angular_dist(ra1,dec1,ra2,dec2)
      if dist<max_dist:
        dists.append(dist)
      else:
        dists.append(5*dist)
    else:
      dists.append(10*max_dist)

  dists = np.array(dists)

  matches = list(zip(1+np.where(dists<=max_dist)[0], 1+idx[dists<=max_dist], dists[dists<=max_dist].tolist()))
  no_matches = (1+np.where((idx==tree.n)|(dists>max_dist))[0]).tolist()
  return matches, no_matches

bss_cat = import_bss()
super_cat = import_super()

# First example in the question
max_dist = 40/3600
matches, no_matches = crossmatch(bss_cat, super_cat, max_dist)
print(matches[:3])
print(no_matches[:3])
print(len(no_matches))

# Second example in the question
max_dist = 5/3600
matches, no_matches = crossmatch(bss_cat, super_cat, max_dist)
print(matches[:3])
print(no_matches[:3])
print(len(no_matches))

[(1, 2, 0.00010988610938711933), (2, 4, 0.0007649845967243825), (3, 5, 0.00020863352870694596)]
[5, 6, 11]
9
[(1, 2, 0.00010988610938711933), (2, 4, 0.0007649845967243825), (3, 5, 0.00020863352870694596)]
[5, 6, 11]
40

Data-driven Astronomy

Right ascension and Declination

Right ascension and Declination#