import re
import numpy as np
from astropy.coordinates import SkyCoord, concatenate
from astroquery.simbad import Simbad
from astroquery.vizier import Vizier
import astropy.units as u
import pandas as pd
from tqdm import tqdm
# ===================== SIMBAD 查询函数 (无需改动) =====================
def query_targets(center, radius_deg, target_type_codes, label):
print(f"正在为 '{label}' 查询 SIMBAD，使用代码: {target_type_codes}")
custom_simbad = Simbad()
custom_simbad.add_votable_fields("otype", "ra(d)", "dec(d)")
custom_simbad.ROW_LIMIT = -1
try:
result = custom_simbad.query_region(center, radius=radius_deg * u.deg)
except Exception as e:
print(f"查询 {label} 时出错：{e}")
return pd.DataFrame()
if result is None: return pd.DataFrame()
df = result.to_pandas()
df.columns = [x.lower() for x in df.columns]
ra_col_name = 'ra_d' if 'ra_d' in df.columns else 'ra'
dec_col_name = 'dec_d' if 'dec_d' in df.columns else 'dec'
if ra_col_name not in df.columns: return pd.DataFrame()
if isinstance(df[ra_col_name].iloc[0], str):
coords = SkyCoord(df[ra_col_name], df[dec_col_name], unit=(u.hourangle, u.deg))
df['ra_deg'], df['dec_deg'] = coords.ra.deg, coords.dec.deg
else:
df['ra_deg'], df['dec_deg'] = df[ra_col_name], df[dec_col_name]
df = df[df['otype'].notna()]
mask_include = df['otype'].apply(lambda x: any(re.search(r'\b' + re.escape(code) + r'\b', x, re.IGNORECASE) for code in target_type_codes))
mask_exclude = df['otype'].str.contains('Dark Cloud|Nebula|MolCld', case=False, na=False)
final_mask = mask_include & ~mask_exclude
df_filtered = df[final_mask].copy()
if df_filtered.empty:
print(f"未找到符合 {label} 类型的天体。")
return pd.DataFrame()
df_filtered['Type'] = label
print(f"查询并过滤 {label} 成功，找到 {len(df_filtered)} 条记录。")
return df_filtered
# ===================== 交叉匹配函数 (已重写以提高鲁棒性) =====================
def crossmatch_vizier(df, catalog_id, col_name, radius_arcsec):
if df.empty:
print(f"输入为空，跳过 {catalog_id} 的匹配。")
return df
print(f"正在与 {catalog_id} 进行交叉匹配...")
# 1. 创建查询坐标对象
coords_query = SkyCoord(ra=df["ra_deg"].values, dec=df["dec_deg"].values, unit=u.deg, frame="icrs")
# 2. 一次性查询整个区域，而不是每个源单独查
vizier = Vizier(columns=["**"], row_limit=-1, timeout=300) # 获取所有列
try:
# 使用查询区域的中心和最大半径进行一次性查询
# 我们用原始的中心和半径，这样更简单
center_coord = SkyCoord(ra=df["ra_deg"].mean() * u.deg, dec=df["dec_deg"].mean() * u.deg)
# 查询半径需要覆盖整个区域，用原始fov_radius加上一点余量
search_radius = (df['ra_deg'].max() - df['ra_deg'].min()) / 2 + 0.1 # 简化的估算
result_tables = vizier.query_region(center_coord, radius=search_radius * u.deg, catalog=catalog_id)
except Exception as e:
print(f"匹配 {catalog_id} 时出错：{e}")
df[col_name] = False
return df
if not result_tables:
print(f"{catalog_id} 中无匹配结果。")
df[col_name] = False
return df
# 3. 将返回的所有表合并成一个大的 astropy Table
# result_tables[0] 通常是包含所有结果的主表
matched_table = result_tables[0]
if len(matched_table) == 0:
print(f"{catalog_id} 中无匹配结果。")
df[col_name] = False
return df
# 4. 从这个大表中创建星表坐标对象
matched_catalog_coords = SkyCoord(ra=matched_table['_RAJ2000'], dec=matched_table['_DEJ2000'], unit=u.deg)
# 5. 使用 astropy 高效的 match_to_catalog_sky 进行匹配
idx, sep2d, _ = coords_query.match_to_catalog_sky(matched_catalog_coords)
# 6. 创建一个布尔掩码，只有当最小距离小于我们的搜索半径时才为True
matched_mask = sep2d < (radius_arcsec * u.arcsec)
df[col_name] = matched_mask
print(f"为 '{col_name}' 匹配到 {matched_mask.sum()} / {len(df)} 个目标。")
return df
# ===================== 主流程 (无需改动) =====================
if __name__ == "__main__":
ra_deg = 287.55
dec_deg = 5.28
fov_radius = 1.414
center_coord = SkyCoord(ra=ra_deg * u.deg, dec=dec_deg * u.deg, frame="icrs")
print(f"查询中心: {center_coord.to_string('hmsdms')}")
print(f"搜索半径: {fov_radius} 度")
target_types = {
"激变变星": ["CV", "Nova", "DN", "DNe", "AM", "IP"],
"X射线双星": ["XRB", "HMXB", "LMXB"],
"大质量双星": ["SB"],
"大质量恒星": ["O", "B"]
}
all_results = []
for label, otypes in tqdm(target_types.items(), desc="查询SIMBAD分类"):
tqdm.write(f"\n--- 开始查询: {label} ---")
df = query_targets(center_coord, fov_radius, otypes, label)
if not df.empty:
all_results.append(df)
if all_results:
final_df = pd.concat(all_results, ignore_index=True)
print(f"\n查询完毕，找到总计 {len(final_df)} 个符合条件的目标。")
print("\n--- 开始多星表交叉匹配 ---")
final_df = crossmatch_vizier(final_df, "I/355/gaiadr3", "in_gaia", radius_arcsec=2)
final_df = crossmatch_vizier(final_df, "II/246/out", "in_2mass", radius_arcsec=5)
final_df = crossmatch_vizier(final_df, "II/336/apass9", "in_apass", radius_arcsec=3)
final_df["only_in_simbad"] = ~(
final_df["in_gaia"] | final_df["in_2mass"] | final_df["in_apass"]
)
save_path = r"D:\谷歌\数据\target_stars_final_robust_v2.csv"
final_df.to_csv(save_path, index=False, encoding="utf-8-sig")
print(f"\n✨ 分析完成！结果已保存至：{save_path}")
print(f"其中，仅在SIMBAD中找到，未匹配到其他星表的目标有 {final_df['only_in_simbad'].sum()} 个。")
else:
print("\n在指定天区内没有找到任何符合条件的天体。")
最新发布
本文介绍了Python pandas库中DataFrame的limit方法，该方法用于获取DataFrame的前n行记录，返回一个新的DataFrame对象。limit不同于take和head，它不是Action操作，意味着它不会立即执行计算，而是保留了操作的延迟特性。
扫一扫
4139
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
