# Week 2 - Crossmatching Catalogues
#### 


```python
import numpy as np
import time
from astropy.coordinates import SkyCoord
from astropy import units as u
```

### Convert from HMS & DMS notation to Decimal degrees


```python
def hms2dec(h,m,s):
  return 15*(h + m/60 + s/(60*60))

def dms2dec(h,m,s):
  return h*(1 + m/(abs(h)*60) + s/(abs(h)*60*60))  

# 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.955000000000002
    -66.08475


#### 

### Haversine Formula for calculating Angular Distance
###### 
\begin{align}
d = 2 \arcsin \sqrt{ \sin^2 \frac{|\delta_1 - \delta_2|}{2} + \cos \delta_1 \cos \delta_2 \sin^2 \frac{|\alpha_1 - \alpha_2|}{2} }
\end{align}


```python
def angular_dist(r1,d1,r2,d2):
  r1, r2, d1, d2 = np.radians(r1), np.radians(r2), np.radians(d1), np.radians(d2)
  a = np.sin(np.abs(d1 - d2)/2)**2
  b = np.cos(d1)*np.cos(d2)*np.sin(np.abs(r1 - r2)/2)**2
  d = 2*np.arcsin(np.sqrt(a + b))
  return np.degrees(d)

# Run your function with the first example in the question.
print(angular_dist(21.07, 0.1, 21.15, 8.2))

# Run your function with the second example in the question
print(angular_dist(10.3, -3, 24.3, -29))
```

    8.100392318146504
    29.208498180546595


#### 

### Reading data from AT20G BSS and SuperCOSMOS catalogues


```python
def import_bss(path):
  data = np.loadtxt(path, usecols=range(1, 7))
  out = []
  for i, row in enumerate(data, 1):
    out.append((i, hms2dec(row[0], row[1], row[2]), dms2dec(row[3], row[4], row[5])))                  
  return out

def import_super(path):
  data = np.loadtxt(path, delimiter=',', skiprows=1, usecols=[0, 1])
  out = []
  for i, row in enumerate(data, 1):
    out.append((i, row[0], row[1]))  
  return out

# Output of the import_bss and import_super functions
bss_cat = import_bss('Data 2/bss_truncated.dat')
super_cat = import_super('Data 2/super_truncated.csv')

print('Object ID  |  Right Ascension°  |  Declination°\n')
print(bss_cat)
print(super_cat)
```

    Object ID  |  Right Ascension°  |  Declination°
    
    [(1, 1.1485416666666666, -47.60530555555555), (2, 2.6496666666666666, -30.463416666666664), (3, 2.7552916666666665, -26.209194444444442)]
    [(1, 1.0583407, -52.9162402), (2, 2.6084425, -41.5005753), (3, 2.7302499, -27.706955)]


#### 

### Finding closest neighbour for a target source (RA°, Dec°) from a catalogue


```python
def find_closest(data, RA1, Dec1):
  ind = 0
  closest = angular_dist(RA1, Dec1, data[0][1], data[0][2])
  for i, row in enumerate(data, 0):
    test = angular_dist(RA1, Dec1, row[1], row[2])
    if test < closest:
      ind = i
      closest = test
  return (data[ind][0], closest)

cat = import_bss('Data 2/bss.dat')
print('ID  |  Angular Distance°\n')

# First example from the question
print(find_closest(cat, 175.3, -32.5))

# Second example in the question
print(find_closest(cat, 32.2, 40.7))

```

    ID  |  Angular Distance°
    
    (156, 3.7670580226469013)
    (26, 57.729135775621295)


#### 

## Crossmatching 2 catalogues within a given distance


```python
def crossmatch(cat1, cat2, dist):
  matches, no_matches = [], []
  for i, row in enumerate(cat1,1):
    test = find_closest(cat2, row[1], row[2])
    if  test[1] < dist:
      matches.append((i, test[0], test[1]))
    else:
      no_matches.append(i)
  return matches, no_matches

bss_cat = import_bss('Data 2/bss (2).dat')
super_cat = import_super('Data 2/super.csv')

# First example in the question
max_dist = 40/3600
matches, no_matches = crossmatch(bss_cat, super_cat, max_dist)
print('1st Object ID  |  2nd Object ID  |  Angular Distance°\n')
print(matches[:3])
print('Unmatched IDs from 1st Catalogue - ', no_matches[:3])
print('No. of Unmatched objects in 1st Catalogue = ', len(no_matches), '\n')

# Second example in the question
max_dist = 5/3600
matches, no_matches = crossmatch(bss_cat, super_cat, max_dist)
print(matches[:3])
print('Unmatched IDs from 1st Catalogue - ', no_matches[:3])
print('No. of Unmatched objects in 1st Catalogue = ', len(no_matches))

```

    1st Object ID  |  2nd Object ID  |  Angular Distance°
    
    [(1, 2, 0.00010988610939332616), (2, 4, 0.0007649845967220993), (3, 5, 0.00020863352870707666)]
    Unmatched IDs from 1st Catalogue -  [5, 6, 11]
    No. of Unmatched objects in 1st Catalogue =  9 
    
    [(1, 2, 0.00010988610939332616), (2, 4, 0.0007649845967220993), (3, 5, 0.00020863352870707666)]
    Unmatched IDs from 1st Catalogue -  [5, 6, 11]
    No. of Unmatched objects in 1st Catalogue =  40


#### 

### Microoptimising the crossmatch


```python
def angular_dist(r1,d1,r2,d2):
  a = np.sin(np.abs(d1 - d2)/2)**2
  b = np.cos(d1)*np.cos(d2)*np.sin(np.abs(r1 - r2)/2)**2
  d = 2*np.arcsin(np.sqrt(a + b))
  return d

def find_closest(data, RA1, Dec1):
  ind = 0
  closest = angular_dist(RA1, Dec1, data[0][0], data[0][1])
  for i, row in enumerate(data, 0):
    test = angular_dist(RA1, Dec1, row[0], row[1])
    if test < closest:
      closest = test
      ind = i
  return (ind, closest)

def crossmatch(cat1, cat2, dist):
  start = time.perf_counter()
  matches, no_matches = [], []
  cat1 = np.radians(cat1)
  cat2 = np.radians(cat2)
  dist = np.radians(dist)
  for i, row in enumerate(cat1,0):
    test = find_closest(cat2, row[0], row[1])
    if  test[1] < dist:
      matches.append((i, test[0], np.degrees(test[1])))
    else:
      no_matches.append(i)
  seconds = time.perf_counter() - start
  return matches, no_matches, seconds


# The example in the question
cat1 = np.array([[180, 30], [45, 10], [300, -45]])
cat2 = np.array([[180, 32], [55, 10], [302, -44]])
matches, no_matches, time_taken = crossmatch(cat1, cat2, 5)
print('1st Object ID  |  2nd Object ID  |  Angular Distance°\n')
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken, '\n')

# A function to create a random catalogue of size n
def create_cat(n):
    ras = np.random.uniform(0, 360, size=(n, 1))
    decs = np.random.uniform(-90, 90, size=(n, 1))
    return np.hstack((ras, decs))

# Test your function on random inputs
np.random.seed(0)
cat1 = create_cat(10)
cat2 = create_cat(20)
matches, no_matches, time_taken = crossmatch(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken)
```

    1st Object ID  |  2nd Object ID  |  Angular Distance°
    
    matches: [(0, 0, 2.0000000000000027), (2, 2, 1.7420109046547023)]
    unmatched: [1]
    time taken: 0.0003623000000061438 
    
    matches: []
    unmatched: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    time taken: 0.005214199999954872


#### 

### Vectorisation using NumPy


```python
def crossmatch_vect(cat1, cat2, max_radius):
  start = time.perf_counter()
  max_radius = np.radians(max_radius)
  
  matches, no_matches = [], []

  # Convert coordinates to radians
  cat1 = np.radians(cat1)
  cat2 = np.radians(cat2)
  ra2s = cat2[:,0]
  dec2s = cat2[:,1]

  for id1, (ra1, dec1) in enumerate(cat1):
    dists = angular_dist(ra1, dec1, ra2s, dec2s)
    min_id = np.argmin(dists)
    min_dist = dists[min_id]
    if min_dist > max_radius:
      no_matches.append(id1)
    else:
      matches.append((id1, min_id, np.degrees(min_dist)))
    
  time_taken = time.perf_counter() - start
  return matches, no_matches, time_taken

# The example in the question
ra1, dec1 = np.radians([180, 30])
cat2 = [[180, 32], [55, 10], [302, -44]]
cat2 = np.radians(cat2)
ra2s, dec2s = cat2[:,0], cat2[:,1]
dists = angular_dist(ra1, dec1, ra2s, dec2s)
print('Angular distance° - ', np.degrees(dists), '\n')

cat1 = np.array([[180, 30], [45, 10], [300, -45]])
cat2 = np.array([[180, 32], [55, 10], [302, -44]])
matches, no_matches, time_taken = crossmatch_vect(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken, '\n')

# Test your function on random inputs
cat1 = create_cat(10)     # Create a random catalogue of size 10
cat2 = create_cat(20)
matches, no_matches, time_taken = crossmatch_vect(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken)
```

    Angular distance° -  [  2.         113.72587199 132.64478705] 
    
    matches: [(0, 0, 2.0000000000000027), (2, 2, 1.7420109046547023)]
    unmatched: [1]
    time taken: 0.00032420000025012996 
    
    matches: []
    unmatched: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    time taken: 0.0007706999999754771


#### 

### Breaking out after maximum match radius : Searching within -90° < δ° < (δ + r)°


```python
def crossmatch(cat1, cat2, max_radius):
  start = time.perf_counter()
  max_radius = np.radians(max_radius)  
  matches, no_matches = [], []

  cat1 = np.radians(cat1)
  cat2 = np.radians(cat2)
  order = np.argsort(cat2[:,1])
  cat2_ordered = cat2[order]
  
  for id1, (ra1, dec1) in enumerate(cat1):
    min_dist = np.inf
    min_id2 = None
    max_dec = dec1 + max_radius
    for id2, (ra2, dec2) in enumerate(cat2_ordered):
      if dec2 > max_dec:
        break
      dist = angular_dist(ra1, dec1, ra2, dec2)
      if dist < min_dist:
        min_id2 = order[id2]
        min_dist = dist
    if min_dist > max_radius:
      no_matches.append(id1)
    else:
      matches.append((id1, min_id2, np.degrees(min_dist)))
    
  time_taken = time.perf_counter() - start
  return matches, no_matches, time_taken


# The example in the question
cat1 = np.array([[180, 30], [45, 10], [300, -45]])
cat2 = np.array([[180, 32], [55, 10], [302, -44]])
matches, no_matches, time_taken = crossmatch(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken, '\n')

# Test your function on random inputs
np.random.seed(0)
cat1 = create_cat(10)
cat2 = create_cat(20)
matches, no_matches, time_taken = crossmatch(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken)
```

    matches: [(0, 0, 2.0000000000000027), (2, 2, 1.7420109046547023)]
    unmatched: [1]
    time taken: 0.015049199999793927 
    
    matches: []
    unmatched: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    time taken: 0.006830000000263681


#### 

### Boxing Match : Searching within (δ - r)° < δ° < (δ + r)°


```python
def crossmatch_boxing(cat1, cat2, max_radius):
  start = time.perf_counter()
  max_radius = np.radians(max_radius)  
  matches, no_matches = [], []

  cat1 = np.radians(cat1)
  cat2 = np.radians(cat2)
  order = np.argsort(cat2[:,1])
  cat2_ordered = cat2[order]
  
  for id1, (ra1, dec1) in enumerate(cat1):
    min_dist = np.inf
    min_id2 = None
    max_dec = dec1 + max_radius
    index = np.searchsorted(cat2_ordered[:,1], dec1 - max_radius, side='left')
    
    for id2, (ra2, dec2) in enumerate(cat2_ordered[index:,:]):
      if dec2 > max_dec:
        break
      dist = angular_dist(ra1, dec1, ra2, dec2)
      if dist < min_dist:
        min_id2 = order[index:][id2]
        min_dist = dist
    if min_dist > max_radius:
      no_matches.append(id1)
    else:
      matches.append((id1, min_id2, np.degrees(min_dist)))
    
  time_taken = time.perf_counter() - start
  return matches, no_matches, time_taken


# The example in the question
cat1 = np.array([[180, 30], [45, 10], [300, -45]])
cat2 = np.array([[180, 32], [55, 10], [302, -44]])
matches, no_matches, time_taken = crossmatch_boxing(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken, '\n')

# Test your function on random inputs
np.random.seed(0)
cat1 = create_cat(10)
cat2 = create_cat(20)
matches, no_matches, time_taken = crossmatch_boxing(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken)
```

    matches: [(0, 0, 2.0000000000000027), (2, 2, 1.7420109046547023)]
    unmatched: [1]
    time taken: 0.044778799999676266 
    
    matches: []
    unmatched: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    time taken: 0.0005314999998518033


#### 

### Crossmatching with k-d Trees


```python
def crossmatch_kd(cat1, cat2, dist):
  start = time.perf_counter()
  sky_cat1 = SkyCoord(cat1*u.degree, frame='icrs')
  sky_cat2 = SkyCoord(cat2*u.degree, frame='icrs')
  closest_ids, closest_dists, closest_dists3d = sky_cat1.match_to_catalog_sky(sky_cat2) 
  matches, no_matches = [], []
  for i, ele in enumerate(closest_dists.value):
    if  ele < dist:
      matches.append((i, closest_ids[i], ele))
    else:
      no_matches.append(i)
  seconds = time.perf_counter() - start
  return matches, no_matches, seconds


# The example in the question
cat1 = np.array([[180, 30], [45, 10], [300, -45]])
cat2 = np.array([[180, 32], [55, 10], [302, -44]])
matches, no_matches, time_taken = crossmatch_kd(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken, '\n')

# Test your function on random inputs
np.random.seed(0)
cat1 = create_cat(10)
cat2 = create_cat(20)
matches, no_matches, time_taken = crossmatch_kd(cat1, cat2, 5)
print('matches:', matches)
print('unmatched:', no_matches)
print('time taken:', time_taken)
```

    matches: [(0, 0, 2.0000000000000036), (2, 2, 1.7420109046547163)]
    unmatched: [1]
    time taken: 2.324950500000341 
    
    matches: []
    unmatched: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    time taken: 0.005834800000229734


a317d319b0aa247782cb2152703ed7babc01bf4e
ipynb
Jupyter Notebook
Week 2 - Crossmatching/Week 2.ipynb
utsav-akhaury/Data-driven-Astronomy
7b09c30054a46f915b1a3e88b0cf59f91a4308f5
MIT
Week 2 - Crossmatching/Week 2.ipynb
utsav-akhaury/Data-driven-Astronomy
7b09c30054a46f915b1a3e88b0cf59f91a4308f5
MIT
Week 2 - Crossmatching/Week 2.ipynb
utsav-akhaury/Data-driven-Astronomy
7b09c30054a46f915b1a3e88b0cf59f91a4308f5
MIT
Qwen/Qwen-72B
1. YES
2. YES