In [1]:
from __future__ import division, print_function # 使程式碼兼容 Python 2
import time # 計時
import numpy as np # 科學運算套件
import pandas as pd # 讀取 csv 檔
In [2]:
table_a = pd.read_csv('files/table_A.csv')
table_b = pd.read_csv('files/table_B.csv')
In [3]:
table_a.head()
Out[3]:
In [4]:
table_b.head()
Out[4]:
In [5]:
print('Table A 有', len(table_a), '個樣本')
print('Table B 有', len(table_b), '個樣本')
In [6]:
ra_a = table_a['ra']
dec_a = table_a['dec']
ra_b = table_b['wise_ra']
dec_b = table_b['wise_dec']
In [7]:
t1_0 = time.time() # 計時器_開始時間
idx1 = []
for i in range(len(table_a[:1])):
dist1 = []
for j in range(len(table_b)):
dist1.append(np.sqrt((ra_a[i] - ra_b[j])**2 + (dec_a[i] - dec_b[j])**2))
if min(dist1) < (1. / 3600.):
idx1.append(dist1.index(min(dist1)))
else:
idx1.append(-9999)
t1_1 = time.time() # 計時器_結束時間
dt1 = t1_1 - t1_0
dt1t = dt1 / (len(table_a[:1]) / len(table_a))
In [8]:
print('配對一個需要花費：', round(dt1, 2), '秒')
print('配對全部需要花費：', round(dt1t / 3600, 2), '小時')
print(idx1)
In [9]:
t15_0 = time.time() # 計時器_開始時間
idx15 = np.zeros_like(ra_a)
dist15 = np.zeros_like(ra_a)
for i in range(len(table_a[:5000])):
ind = np.where(np.sqrt((ra_b - ra_a[i])**2 + (dec_b - dec_a[i])**2) < (1. / 3600.))[0]
if ind.size == 0:
idx15[i] = -9999
dist15[i] = -9999
else:
idx15[i] = ind[0]
dist15[i] = np.sqrt((ra_b[ind[0]] - ra_a[i])**2 + (dec_b[ind[0]] - dec_a[i])**2)
t15_1 = time.time() # 計時器_結束時間
dt15 = t15_1 - t15_0
dt15t = dt15 / (len(table_a[:5000]) / len(table_a))
In [10]:
print('配對5000個需要花費：', round(dt15, 2), '秒')
print('配對全部需要花費：', round(dt15t, 2), '秒')
print(idx15[:20])
In [11]:
from astropy.coordinates import SkyCoord # astropy的座標套件
from astropy import units as u # astropy的單位套件
In [12]:
t2_0 = time.time() # 計時器_開始時間
# 座標輸入
ca = SkyCoord(ra_a*u.deg, dec_a*u.deg)
cb = SkyCoord(ra_b*u.deg, dec_b*u.deg)
# 配對並計算相差距離
idx2, d2d, d3d = ca.match_to_catalog_sky(cb)
matches = cb[idx2]
dra, ddec = ca.spherical_offsets_to(matches)
# 設定配對
idx2_new = idx2[d2d < 5 * u.arcsec]
t2_1 = time.time() # 計時器_結束時間
dt2 = t2_1 - t2_0
In [13]:
print('配對全部需要花費：', round(dt2, 2), '秒')
print(idx2_new[:20])
In [14]:
from astroML.crossmatch import crossmatch_angular
In [15]:
t3_0 = time.time() # 計時器_開始時間
# get imaging data
imX = np.empty((len(table_a), 2), dtype=np.float64)
imX[:, 0] = ra_a
imX[:, 1] = dec_a
# get standard stars
stX = np.empty((len(table_b), 2), dtype=np.float64)
stX[:, 0] = ra_b
stX[:, 1] = dec_b
# crossmatch catalogs
max_radius = 1. / 3600 # 1 arcsec
dist3, idx3 = crossmatch_angular(imX, stX, max_radius)
match = ~np.isinf(dist3)
# 將沒有 match 到的設為 -9999
idx3[~match] = -9999
# 配對距離 (單位：角秒)
dist_match = dist3[match]
dist_match *= 3600
t3_1 = time.time() # 計時器_結束時間
dt3 = t3_1 - t3_0
In [16]:
print('配對全部需要花費：', round(dt3, 2), '秒')
print(idx3[:20])
In [17]:
print('方法一：簡單交叉配對 (Naive cross-match)\n\n', '\t需耗時：', round(dt1t / 3600, 2), '小時\n')
print('方法1.5：改用 numpy array\n\n', '\t需耗時：', round(dt15t, 2), '秒\n')
print('方法二：Astropy 交叉配對\n\n', '\t需耗時：', round(dt2, 2), '秒\n')
print('方法三：Astroml 交叉配對\n\n', '\t需耗時：', round(dt3, 2), '秒\n')