如何用Python基于笛卡尔坐标交叉匹配两个天文目录DataFrame？

优先推荐使用Astropy完成这项工作——它专门针对天文天球坐标的球面几何匹配做了优化，比纯numpy/pandas的平面匹配精度更高，完全适配巡天星表的匹配需求，具体步骤如下：

1. 导入依赖库 #

import pandas as pd
from astropy.coordinates import SkyCoord
from astropy import units as u

2. 加载星表数据 #

假设你的两个星表已经是pandas DataFrame格式（比如df1和df2），核心列包含ra（赤经）和dec（赤纬）：

# 替换为你的数据读取逻辑，比如从CSV/FITS加载
df1 = pd.read_csv('catalog1.csv')
df2 = pd.read_csv('catalog2.csv')

3. 转换为天球坐标对象 #

把DataFrame中的坐标转换为Astropy的SkyCoord类型，支持直接处理球面几何：

# 情况1：坐标单位为度（Decimal Degrees）
coord1 = SkyCoord(ra=df1['ra']*u.degree, dec=df1['dec']*u.degree)
coord2 = SkyCoord(ra=df2['ra']*u.degree, dec=df2['dec']*u.degree)

# 情况2：坐标为时分秒/度分秒格式（比如ra是'12h34m56s'，dec是'+45d30m15s'）
# coord1 = SkyCoord(ra=df1['ra'], dec=df1['dec'], unit=(u.hourangle, u.degree))

4. 执行交叉匹配 #

使用match_to_catalog_sky方法完成球面匹配，设置合理的匹配半径（比如2角秒，根据你的巡天观测精度调整）：

# 匹配得到：df2中与df1每个源最接近的索引、角距离、3D距离
idx, d2d, d3d = coord1.match_to_catalog_sky(coord2)

# 筛选出符合距离阈值的匹配项
match_threshold = 2 * u.arcsec  # 可根据需求修改
match_mask = d2d < match_threshold

# 提取匹配上的两行数据并重置索引
matched_df1 = df1[match_mask].reset_index(drop=True)
matched_df2 = df2[idx[match_mask]].reset_index(drop=True)

# 合并成最终的交叉匹配星表，给第二个星表的列加前缀区分
final_catalog = pd.concat([matched_df1, matched_df2.add_prefix('cat2_')], axis=1)
# 新增匹配距离列，单位转换为角秒
final_catalog['match_sep_arcsec'] = d2d[match_mask].to_value(u.arcsec)

5. 处理一对多/多对一匹配（可选） #

如果出现一个源匹配到多个目标的情况，可以筛选每个源的最近匹配：

# 对df2的每个匹配目标，找到df1中距离最近的源
min_dist_idx = d2d.groupby(idx).idxmin()
unique_matched_df1 = df1.iloc[min_dist_idx].reset_index(drop=True)
unique_matched_df2 = df2.iloc[idx[min_dist_idx]].reset_index(drop=True)

备选方案：纯Pandas+SciPy粗略匹配（不推荐高精度场景） #

如果不想引入Astropy，可以将天球坐标转换为笛卡尔坐标后用KDTree匹配，但球面转平面会存在投影误差，仅适合低精度需求：

import numpy as np
from scipy.spatial import cKDTree

# 天球坐标转笛卡尔坐标函数
def sky_to_cartesian(ra, dec):
    ra_rad = np.radians(ra)
    dec_rad = np.radians(dec)
    x = np.cos(dec_rad) * np.cos(ra_rad)
    y = np.cos(dec_rad) * np.sin(ra_rad)
    z = np.sin(dec_rad)
    return np.array([x, y, z]).T

# 生成坐标数组
cart1 = sky_to_cartesian(df1['ra'], df1['dec'])
cart2 = sky_to_cartesian(df2['ra'], df2['dec'])

# KDTree匹配
tree = cKDTree(cart2)
distances, indices = tree.query(cart1, k=1)

# 转换为角距离（弧度转角秒）
angle_dist = 2 * np.arcsin(distances / 2) * (180/np.pi)*3600
match_mask = angle_dist < 2

# 合并生成最终星表
final_catalog = pd.concat([df1[match_mask], df2.iloc[indices[match_mask]].add_prefix('cat2_')], axis=1)
final_catalog['match_sep_arcsec'] = angle_dist[match_mask]

关键提示：优先选择Astropy方案，它完全适配天球坐标的球面特性，不会因投影畸变导致匹配错误，是天文星表匹配的标准工具。

内容的提问来源于stack exchange，提问作者NeStack