Skip to content

Instantly share code, notes, and snippets.

@ArnulfoPerez
Forked from lucassouzavieira/cross_matching.py
Created December 13, 2020 21:10
Show Gist options
  • Select an option

  • Save ArnulfoPerez/1ab397578e827ae37559a213ab3b1988 to your computer and use it in GitHub Desktop.

Select an option

Save ArnulfoPerez/1ab397578e827ae37559a213ab3b1988 to your computer and use it in GitHub Desktop.
import numpy as np
# Angular distance
def angular_distance(r1, d1, r2, d2):
r1 = np.radians(r1)
r2 = np.radians(r2)
d1 = np.radians(d1)
d2 = 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
angular_distance = 2 * np.arcsin(np.sqrt(a + b))
return np.degrees(angular_distance)
# HMS2dec
def hms2dec(hours, minutes, seconds):
angle = 15 * (hours + minutes / 60 + seconds / (60 * 60))
return angle
# DMS2dec
def dms2dec(degrees, arcmin, arcsec):
negation = 1
if degrees < 0:
degrees = degrees * -1
negation = -1
angle = negation * (degrees + arcmin / 60 + arcsec/(60*60))
return angle
# Import BSS
def import_bss():
data = []
bss = np.loadtxt("bss.dat", usecols=range(1,7))
for index, line in enumerate(bss):
data.append((index + 1, hms2dec(line[0], line[1], line[2]), dms2dec(line[3], line[4], line[5])))
return data
# Import Super
def import_super():
data = []
catalogue = np.loadtxt('super.csv', delimiter=',', skiprows=1,usecols=[0,1])
for index, line in enumerate(catalogue):
data.append((index + 1, line[0], line[1]))
return data
# Find the closest match
def find_closest(dataset, right_ascension, declination):
min_distance, min_id = np.inf, 0
for id1, r1, d1 in dataset:
distance = angular_distance(r1, d1, right_ascension, declination)
if min_distance > distance:
min_id = id1
min_distance = distance
return min_id, min_distance
# Cross-match function
def crossmatch(dataset1, dataset2, max_radius):
matches = []
no_matches = []
for id1, ra1, dec1 in dataset1:
closest_dist = np.inf
closest_id2 = None
for id2, ra2, dec2 in dataset2:
dist = angular_distance(ra1, dec1, ra2, dec2)
if dist < closest_dist:
closest_id2 = id2
closest_dist = dist
# Ignore match if it's outside the maximum radius
if closest_dist > max_radius:
no_matches.append(id1)
else:
matches.append((id1, closest_id2, closest_dist))
return matches, no_matches
if __name__ == '__main__':
cat = import_bss()
# 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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment