In [1]:
from __future__ import division, print_function # 使程式碼兼容 Python 2
import time # 計時
import numpy as np # 科學運算套件
import pandas as pd # 讀取 csv 檔

載入 tables


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]:
ra dec z specObjID objid u_mag g_mag r_mag i_mag z_mag
0 187.27970 -1.193372 0.032031 3.250000e+17 1.240000e+18 19.24800 18.01970 17.60396 17.41241 17.15981
1 210.76184 -1.094643 0.032493 3.400000e+17 1.240000e+18 18.59815 17.87038 17.66903 17.52236 17.59061
2 185.74969 -0.293052 0.030245 3.250000e+17 1.240000e+18 19.02032 18.05203 17.60972 17.43984 17.47915
3 193.52240 0.175311 0.022748 3.290000e+17 1.240000e+18 18.30348 16.88793 16.19267 15.73816 15.40265
4 222.84758 0.007093 0.023165 3.480000e+17 1.240000e+18 17.55991 15.81113 15.03432 14.67037 14.45513

In [4]:
table_b.head()


Out[4]:
wise_ra wise_dec sdss_objid w1 w1err w1snr w2 w2err w2snr w3 w3err w3snr w4 w4err w4snr
0 0.697183 29.359842 1237663307454874339 14.330 0.030 36.7 14.054 0.046 23.6 12.702 9999.000 -1.5 9.230 9999.0 -1.1
1 359.329741 32.799931 1237663306916430563 17.363 0.190 5.7 17.097 9999.000 -0.2 12.429 9999.000 1.0 9.065 9999.0 0.5
2 1.379313 24.340139 1237663306920231248 17.354 0.177 6.1 16.286 0.314 3.5 12.748 9999.000 -0.6 8.721 9999.0 1.0
3 2.291314 24.765421 1237663307993842580 16.353 0.074 14.7 16.233 0.281 3.9 12.498 9999.000 -0.6 8.351 9999.0 1.6
4 2.475097 26.254814 1237663308530122881 13.670 0.027 40.8 13.599 0.038 28.6 11.877 0.246 4.4 9.057 9999.0 0.4

In [5]:
print('Table A 有', len(table_a), '個樣本')
print('Table B 有', len(table_b), '個樣本')


Table A 有 59574 個樣本
Table B 有 129350 個樣本

In [6]:
ra_a = table_a['ra']
dec_a = table_a['dec']

ra_b = table_b['wise_ra']
dec_b = table_b['wise_dec']

方法一:簡單交叉配對 (Naive cross-match)

使用兩層 for 迴圈,並設定 5 角秒的配對距離

  • 距離 $= \sqrt{(RA_{ai}-RA_{bi})^{2}+(Dec_{ai}-Dec_{bi})^{2}}$
  • idx 為 table_a 對應到 table_b 的 index
  • 例如 table_a 的第1列配對到 table_b 的第53列,idx 就是 53

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)


配對一個需要花費: 13.9 秒
配對全部需要花費: 230.02 小時
[62698]

方法1.5:改用 Numpy array


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])


配對5000個需要花費: 20.05 秒
配對全部需要花費: 238.85 秒
[  62698.   -9999.   61225.   -9999.   81327.   87556.   89451.   -9999.
   -9999.  111927.  122692.   -9999.   58221.   70848.  105645.   -9999.
   59439.   69199.   73789.  105999.]

方法二:Astropy 交叉配對


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])


配對全部需要花費: 17.12 秒
[ 62698 101885  61225  67396  81327  87556  89451 121967 124333 111927
 122692  58221  70848 105645  50528  59439  69199  73789 105999  86699]

方法三:Astroml 交叉配對


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])


配對全部需要花費: 0.21 秒
[ 62698  -9999  61225  -9999  81327  87556  89451  -9999  -9999 111927
 122692  -9999  58221  70848 105645  -9999  59439  69199  73789 105999]

運算時間


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')


方法一:簡單交叉配對 (Naive cross-match)

 	需耗時: 230.02 小時

方法1.5:改用 numpy array

 	需耗時: 238.85 秒

方法二:Astropy 交叉配對

 	需耗時: 17.12 秒

方法三:Astroml 交叉配對

 	需耗時: 0.21 秒

方法說明

方法一:使用 Python 原生功能,表現極差

方法1.5:改用 Numpy array 重新撰寫,改善 99.9 %

方法二:使用 Astropy 套件,但 SkyCoord 這個步驟需要花費比較多時間,相較方法1.5,改善 92.8 %

方法三:使用 Astroml 套件,配對方式與方法二同樣為 KDtree,但不須執行 SkyCoord 轉換,相較方法二,改善 98.8 %,在1秒內即可完成