如何在Astropy中围绕原点正确旋转Quadrangle？

要在天图上正确旋转Quadrangle（模拟狭缝），核心是区分天球空间旋转和像素空间旋转，前者能保证狭缝在天球上的几何正确性，避免投影畸变。以下是两种可靠的实现方法：

方法一：天球空间定义旋转后顶点（推荐，无畸变） #

这种方法直接在天球坐标系中计算旋转后的狭缝顶点，再投影到像素空间，适用于所有WCS投影类型，不会产生畸变。

import matplotlib.pyplot as plt
import numpy as np
from astropy.wcs import WCS
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
from astropy import units as u
from astropy.coordinates import SkyCoord
from matplotlib.patches import Polygon

# 加载天图数据
filename = get_pkg_data_filename('galactic_center/gc_msx_e.fits')
hdu = fits.open(filename)[0]
wcs = WCS(hdu.header)

# 创建WCS投影轴
ax = plt.subplot(projection=wcs)
ax.imshow(hdu.data, vmin=-2.e-5, vmax=2.e-4, origin='lower')

# 定义狭缝参数
center_sky = SkyCoord(ra=266.15*u.deg, dec=-28.825*u.deg, frame='fk5')  # 原Quadrangle的中心
slit_width = 0.3*u.deg       # 赤经方向宽度
slit_height = 0.15*u.deg     # 赤纬方向高度
rot_angle = 30*u.deg         # 逆时针旋转角度

# 1. 在局部笛卡尔坐标系（以狭缝中心为原点）定义原始顶点
local_verts = np.array([
    [-slit_width.value/2, -slit_height.value/2],
    [slit_width.value/2, -slit_height.value/2],
    [slit_width.value/2, slit_height.value/2],
    [-slit_width.value/2, slit_height.value/2]
])

# 2. 应用旋转矩阵
theta = np.deg2rad(rot_angle.value)
rot_matrix = np.array([
    [np.cos(theta), -np.sin(theta)],
    [np.sin(theta), np.cos(theta)]
])
rotated_local_verts = local_verts @ rot_matrix.T

# 3. 转换为天球坐标（小角度近似，狭缝尺寸足够小时精度足够）
ra_offsets = rotated_local_verts[:, 0] * u.deg / np.cos(center_sky.dec.rad)
dec_offsets = rotated_local_verts[:, 1] * u.deg
rotated_verts_sky = SkyCoord(
    ra=center_sky.ra + ra_offsets,
    dec=center_sky.dec + dec_offsets,
    frame='fk5'
)

# 4. 转换为像素坐标并绘制Polygon
rotated_verts_pix = wcs.world_to_pixel(rotated_verts_sky)
slit_poly = Polygon(rotated_verts_pix, edgecolor='blue', facecolor='none', lw=2, label='Rotated Slit')
ax.add_patch(slit_poly)

# 绘制原始Quadrangle作为对比
from astropy.visualization.wcsaxes import Quadrangle
original_q = Quadrangle((266.0, -28.9)*u.deg, 0.3*u.deg, 0.15*u.deg,
                        edgecolor='green', facecolor='none', label='Original Quadrangle')
ax.add_patch(original_q)

ax.set_xlabel('Galactic Longitude')
ax.set_ylabel('Galactic Latitude')
plt.legend()
plt.show()

方法二：像素空间旋转（仅适用于线性投影） #

如果你的WCS投影是线性的（如正射投影），可以直接在像素空间围绕狭缝中心旋转，但非线性投影会导致畸变。

import matplotlib.pyplot as plt
import matplotlib.transforms as transforms
import numpy as np
from astropy.wcs import WCS
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
from astropy import units as u
from astropy.visualization.wcsaxes import Quadrangle
from astropy.coordinates import SkyCoord

filename = get_pkg_data_filename('galactic_center/gc_msx_e.fits')
hdu = fits.open(filename)[0]
wcs = WCS(hdu.header)

ax = plt.subplot(projection=wcs)
ax.imshow(hdu.data, vmin=-2.e-5, vmax=2.e-4, origin='lower')

# 定义Quadrangle参数
corner = (266.0, -28.9)*u.deg
width = 0.3*u.deg
height = 0.15*u.deg
rot_angle = 30

# 计算狭缝中心的像素坐标
center_sky = SkyCoord(ra=corner[0]+width/2, dec=corner[1]+height/2, frame='fk5')
center_pix = wcs.world_to_pixel(center_sky)

# 创建围绕中心像素的旋转变换
rot_transform = transforms.Affine2D().rotate_deg_around(*center_pix, rot_angle)
# 获取WCS世界坐标到像素的变换
wcs_transform = ax.get_transform('fk5')

# 组合变换：先转像素再旋转
q = Quadrangle(corner, width, height, edgecolor='purple', facecolor='none', lw=2, label='Pixel-Rotated Quadrangle')
q.set_transform(wcs_transform + rot_transform)
ax.add_patch(q)

ax.set_xlabel('Galactic Longitude')
ax.set_ylabel('Galactic Latitude')
plt.legend()
plt.show()

你之前代码的问题分析 #

1. 旋转原点计算错误：你用r.get_patch_transform().transform((0.5, 0.5))获取的是patch局部坐标（0-1范围）的像素位置，未结合WCS变换，导致原点坐标错误。正确做法是先计算狭缝中心的世界坐标，再转换为像素坐标。
2. 变换顺序错误：td+tr的顺序是先转像素再旋转，对于非线性投影，这种操作会使狭缝形状偏离天球上的实际旋转效果，产生畸变。

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