如何用astropy和regions在指定天区绘制散点图并添加圆形

以下是使用astropy和regions库实现需求的完整代码，支持批量处理多个数据点，同时考虑天区球面曲率：

步骤1：导入依赖库 #

import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord, FK5
from regions import CircleSkyRegion
from astropy.wcs import WCS
from astropy import units as u

步骤2：定义数据与区域参数 #

# 数据点赤道坐标（单位：度）
ra = np.array([248.3, 249.8, 250.2, 250.5, 250.5])
dec = np.array([68.8, 67.7, 65.7, 72.2, 63.3])
# 批量转换为天球坐标对象（支持数组输入，自动处理多个点）
data_points = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame=FK5)

# 圆形区域参数
circle_center = SkyCoord(ra=270*u.degree, dec=66*u.degree, frame=FK5)
circle_radius = 2*u.degree  # 可根据需求调整半径
circle_region = CircleSkyRegion(center=circle_center, radius=circle_radius)

步骤3：创建WCS投影（处理天区曲率） #

为了正确展示天区球面的曲率，我们需要创建一个世界坐标系（WCS），这里选择常用的切平面投影（TAN）：

# 构建WCS参数，以数据点+圆形的覆盖区域为中心设置视场
# 先计算数据点的中心坐标作为视场中心
field_center = SkyCoord(ra=np.mean(ra)*u.degree, dec=np.mean(dec)*u.degree, frame=FK5)
wcs = WCS(naxis=2)
wcs.wcs.crval = [field_center.ra.degree, field_center.dec.degree]  # 视场中心坐标
wcs.wcs.crpix = [150, 150]  # 图像中心像素位置
wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]  # 投影类型
wcs.wcs.cdelt = [-0.1, 0.1]  # 像素对应的角大小（单位：度/像素，负号表示RA从左到右递减）

步骤4：绘制图像 #

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(projection=wcs)

# 绘制数据散点图
ax.scatter(data_points.ra.degree, data_points.dec.degree, 
           transform=ax.get_transform('fk5'),  # 使用FK5赤道坐标系转换
           s=50, c='red', marker='o', label='数据点')

# 绘制圆形区域
circle_pix = circle_region.to_pixel(ax.wcs)  # 转换为像素坐标
circle_pix.plot(ax=ax, facecolor='none', edgecolor='blue', linewidth=2, label='圆形区域')

# 设置坐标轴与图例
ax.set_xlabel('Right Ascension (deg)')
ax.set_ylabel('Declination (deg)')
ax.grid(True, linestyle='--', alpha=0.7)
ax.legend()

plt.show()

关键说明 #

• 多个数据点的处理：SkyCoord直接支持传入数组类型的RA/DEC，无需循环单个创建，批量转换后即可直接用于绘图。
• 天区曲率处理：通过WCS投影将球面坐标转换为平面像素坐标，CircleSkyRegion会自动根据WCS计算符合天区曲率的圆形边界，而非简单的平面圆。
• 视场调整：如果数据点或圆形超出当前视场，可以调整wcs.wcs.cdelt（像素角大小）或crval（视场中心）来适配。

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