问题描述
我有一系列.fit的CCD望远镜观测图，我使用 DAOStarFinder in photutils进行测光。我想找出相同恒星的数据并绘制其光曲线。从一个图形获得的stars表看起来像：
id xcentroid ycentroid sharpness roundness1 roundness2 npix sky peak flux mag JD-HELIO
1 7.2507764 7.1405824 0.25877792 -0.27655266 -0.014479976 361 0 422 12.624582 -2.7530425 2458718.4
2 2740.913 7.1539345 0.30025149 0.25451119 0.018093909 361 0 458 13.467545 -2.8232211 2458718.4
3 2515.323 43.807826 0.44318308 -0.69941143 0.2041087 361 0 626 5.5224866 -1.8553367 2458718.4
4 1284.7828 53.552295 0.980886 -0.1667419 -0.59773071 361 0 480 1.0122505 -0.013219984 2458718.4
5 1249.1764 149.95127 0.49822284 0.27794038 0.72191121 361 0 463 1.0285458 -0.030559111 2458718.4
6 2376.8957 202.02747 0.82856798 -0.42128731 0.30015708 361 0 653 4.7224348 -1.6854149 2458718.4
7 59.982356 215.6136 0.54354939 -0.65521898 -0.00056722803 361 0 497 1.6563912 -0.5479073 2458718.4
8 1661.3238 220.14951 0.42552567 -0.6035932 -0.089775238 361 0 686 7.9859189 -2.2558122 2458718.4
9 1601.4266 229.47563 0.34847327 -0.60959051 0.10700081 361 0 2380 68.697859 -4.592358 2458718.4
10 2559.8217 236.60404 0.2231161 -0.013164036 0.54360116 361 0 490 2.1546284 -0.83343095 2458718.4
... ... ... ... ... ... ... ... ...
Length = 104 rows
第二个数字的表格如下：
id xcentroid ycentroid sharpness roundness1 roundness2 npix sky peak flux mag JD-HELIO
1 7.157499 7.0229365 0.29484857 -0.25416309 -0.018748812 361 0 458 12.940807 -2.7799034 2458718.4
2 2740.8166 7.2593662 0.27311183 0.24388127 0.029043824 361 0 448 12.755378 -2.7642333 2458718.4
3 2495.8085 45.142195 0.64182883 -0.48598172 0.70072363 361 0 649 4.7000496 -1.6802561 2458718.4
4 2679.6403 62.327031 0.80391292 -0.62284532 0.13309338 361 0 524 1.6688956 -0.55607294 2458718.4
5 2467.7183 101.15541 0.64006186 -0.14660082 -0.6205186 361 0 500 1.465029 -0.41461552 2458718.4
6 1094.7268 107.11982 0.59849809 -0.62182068 0.27247315 361 0 476 1.051774 -0.054806086 2458718.4
7 28.506992 145.07421 0.64028434 0.2210654 0.012627233 361 0 474 1.0149705 -0.016133536 2458718.4
8 1876.9044 186.59285 0.82846693 -0.37091888 0.14646202 361 0 515 1.3318741 -0.31115797 2458718.4
9 2358.642 202.09716 0.51573138 -0.31333346 0.88458938 361 0 588 3.6688261 -1.4113178 2458718.4
10 1641.9716 221.59273 0.62435574 -0.56225884 0.25546345 361 0 715 6.4287534 -2.0203169 2458718.4
... ... ... ... ... ... ... ... ...
Length = 84 rows
我对多个图形进行了相同的处理，有趣的是，由于同一颗恒星的xcentroid和ycentroid有一些偏移，如何从不同图形中识别同一颗恒星？ >
找到解决方案
我遇到了同样的问题，并在comaring the CCD stars observation figures中找到了解决方案，其基本思想是找到两组点的三角形的最佳匹配。
然后我用astroalign package计算变换矩阵，并对准所有点。感谢上帝，效果很好。
import itertools
import numpy as np
import matplotlib.pyplot as plt
import astroalign as aa
def getTriangles(set_X,X_combs):
"""
Inefficient way of obtaining the lengths of each triangle's side.
Normalized so that the minimum length is 1.
"""
triang = []
for p0,p1,p2 in X_combs:
d1 = np.sqrt((set_X[p0][0] - set_X[p1][0]) ** 2 +
(set_X[p0][1] - set_X[p1][1]) ** 2)
d2 = np.sqrt((set_X[p0][0] - set_X[p2][0]) ** 2 +
(set_X[p0][1] - set_X[p2][1]) ** 2)
d3 = np.sqrt((set_X[p1][0] - set_X[p2][0]) ** 2 +
(set_X[p1][1] - set_X[p2][1]) ** 2)
d_min = min(d1,d2,d3)
d_unsort = [d1 / d_min,d2 / d_min,d3 / d_min]
triang.append(sorted(d_unsort))
return triang
def sumTriangles(ref_triang,in_triang):
"""
For each normalized triangle in ref,compare with each normalized triangle
in B. find the differences between their sides,sum their absolute values,and select the two triangles with the smallest sum of absolute differences.
"""
tr_sum,tr_idx = [],[]
for i,ref_tr in enumerate(ref_triang):
for j,in_tr in enumerate(in_triang):
# Absolute value of lengths differences.
tr_diff = abs(np.array(ref_tr) - np.array(in_tr))
# Sum the differences
tr_sum.append(sum(tr_diff))
tr_idx.append([i,j])
# Index of the triangles in ref and in with the smallest sum of absolute
# length differences.
tr_idx_min = tr_idx[tr_sum.index(min(tr_sum))]
ref_idx,in_idx = tr_idx_min[0],tr_idx_min[1]
print("Smallest difference: {}".format(min(tr_sum)))
return ref_idx,in_idx
set_ref = np.array([[2511.268821,44.864124],[2374.085032,201.922566],[1619.282942,216.089335],[1655.866502,221.127787],[ 804.171659,2133.549517],])
set_in = np.array([[1992.438563,63.727282],[2285.793346,255.402548],[1568.915358,279.144544],[1509.720134,289.434629],[1914.255205,349.477788],[2370.786382,496.026836],[ 482.702882,508.685952],[2089.691026,523.18825 ],[ 216.827439,561.807396],[ 614.874621,2007.304727],[1286.639124,2155.264827],[ 729.566116,2190.982364]])
# All possible triangles.
ref_combs = list(itertools.combinations(range(len(set_ref)),3))
in_combs = list(itertools.combinations(range(len(set_in)),3))
# Obtain normalized triangles.
ref_triang,in_triang = getTriangles(set_ref,ref_combs),getTriangles(set_in,in_combs)
# Index of the ref and in triangles with the smallest difference.
ref_idx,in_idx = sumTriangles(ref_triang,in_triang)
# Indexes of points in ref and in of the best match triangles.
ref_idx_pts,in_idx_pts = ref_combs[ref_idx],in_combs[in_idx]
print ('triangle ref %s matches triangle in %s' % (ref_idx_pts,in_idx_pts))
print ("ref:",[set_ref[_] for _ in ref_idx_pts])
print ("input:",[set_in[_] for _ in in_idx_pts])
ref_pts = np.array([set_ref[_] for _ in ref_idx_pts])
in_pts = np.array([set_in[_] for _ in in_idx_pts])
transf,(in_list,ref_list) = aa.find_transform(in_pts,ref_pts)
transf_in = transf(set_in)
print(f'transformation matrix: {transf}')
plt.scatter(set_ref[:,0],set_ref[:,1],s=100,marker='.',c='r',label='Reference')
plt.scatter(set_in[:,set_in[:,c='b',label='Input')
plt.scatter(transf_in[:,transf_in[:,marker='+',label='Input Aligned')
plt.plot(ref_pts[:,ref_pts[:,c='r')
plt.plot(in_pts[:,in_pts[:,c='b')
plt.legend()
plt.tight_layout()
plt.savefig( 'align_coordinates.png',format = 'png')
plt.show()
解决方法
当前，我使用一种简单的方法来处理它，如果坐标处于容许范围内，则将其附加。
if (xcentroid-offset < x < xcentroid +offset) and (ycentroid-offset < y < ycentroid+offset):
onestar.append(data[j,:])
实用的整个脚本如下：
#===================================================================
#
# STARS SORTING
#
# This script is used to sorting stars from the table file of each figure.
#
# WARNING:
#
# VERSION: 24 Sep 2020
# AUTHOER: QIANG CHEN chen@camk.edu.pl
#
# PYTHON ENVIRONMENT REQUIREMENT:
# - pip install astropy
# - pip install --no-deps ccdproc
# - pip install photutils
# alternatively can use (CAMK):
# source /home/chen/python-chen/chen3.6/bin/activate
#
# REFERENCE:
# - COMPARISION STARS https://nbviewer.jupyter.org/gist/dokeeffe/416e214f134d39db696c7fdac261964b
#===================================================================
import os
import glob
import numpy as np
PATH = '/work/chuck/chen/obs'
folder = '20190822/reduced'
name = 'ap'
name = 'stars'
offset = 50
filters = ['V','U','B','I']
#filters = ['U']
for filt in filters:
print('SORTING STARS BY EVERY FILTER TYPE:',filt)
files = glob.glob(f'{PATH}/{folder}/{name}_light_*{filt}*.dat')
print(' reading data,setting the most stars file as a reference frame')
line_ref = 0
datas = []
data_ref = np.array([])
for fil in files:
f = open(fil,'r')
next(f) # skip header line
data = np.array([[float(data) for data in line.split()] for line in f.readlines()])
f.close()
datas.append(data)
row,col = data.shape
if (row >= line_ref):
line_ref = row
data_ref = data
if (data_ref.size):
pass
else:
print(' NO RESULTS')
continue
print(' pick up each star in the reference frame,compare x y within all the files')
ii = 0
for i in range(row):
onestar = []
print(' row',i+1,'/',row)
xcentroid = data_ref[i,1]
ycentroid = data_ref[i,2]
for data in datas:
row,col = data.shape
for j in range(row):
x = data[j,1]
y = data[j,2]
if (xcentroid-offset < x < xcentroid +offset) and (ycentroid-offset < y < ycentroid+offset):
onestar.append(data[j,:])
if(len(onestar)<10):
continue
else:
ii+=1
onestar = np.array(onestar)
onestar = onestar[onestar[:,-1].argsort()]
np.savetxt(f'{PATH}/{folder}/sorted_{name}_{filt}{ii}.dat',onestar,fmt='%f ',newline='\n')
设置时间 控制面板
错误1：Request method ‘DELETE‘ not supported 错误还原：...