0

我有一系列.fit的CCD望远镜观测图，我使用photoutils中的DAOStarFinder进行测光。我想找出相同恒星的数据并绘制它们的光变曲线。从一张图中得到的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有一些偏移量，我如何从不同的图形中识别出同一颗星星呢？

找到的解决方案

我遇到了同样的问题，并在 CCD 星星观测图的比较中找到了解决方案，基本思想是找到两组点的三角形的最佳匹配。

然后我使用astroalign包计算变换矩阵，并对齐所有点。感谢主，效果不错。

在此处输入图像描述

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[:,0],set_in[:,1], s=100,marker='.', c='b',label='Input')
plt.scatter(transf_in[:,0],transf_in[:,1], s=100,marker='+', c='b',label='Input Aligned')
plt.plot(ref_pts[:,0],ref_pts[:,1], c='r')
plt.plot(in_pts[:,0],in_pts[:,1], c='b')
plt.legend()
plt.tight_layout()
plt.savefig( 'align_coordinates.png', format = 'png')
plt.show()
4

1 回答 1

0

目前我使用一种简单的方法来处理它，如果坐标在公差范围内，则附加它。

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')
于 2020-09-26T19:52:20.127 回答