如何用imshow实现RA DEC坐标FITS图像转银道坐标系无冗余像素

问题背景 #

我尝试将带有RA/DEC世界坐标（WCS）的FITS图像绘制到银道坐标系中，代码如下：

import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS

fig = plt.figure()
path_galactic = "w40_c18o.fits"
wcs_galactic = WCS(fits.open(path_galactic)[0].header).celestial[0:2] # 参考银道坐标系

path_sofia = 'F0248_FI_IFS_8700081_RED_WXY_100061-100115.fits'
wcs_sofia = WCS(fits.open(path_sofia)[1].header).celestial[0:2] # RA/DEC坐标系的目标图像

a = fig.add_axes([0, 0, 1, 1], projection=wcs_galactic) # 设置绘图投影为银道系
a.set_aspect("equal")
a.imshow(fits.open(path_sofia)[1].data[22],transform=a.get_transform(wcs_sofia)) # 转换坐标系后绘制图像

使用imshow时，图像坐标系转换正确，但出现大量多余无效像素；换成contourf则能得到无冗余的预期效果，需要让imshow实现同等显示效果。

核心原因 #

imshow默认会将图像的矩形像素网格强行填充到整个坐标轴范围，当转换到非笛卡尔坐标系（如银道坐标系）时，原矩形区域外的像素会被拉伸填充，产生冗余。而contourf仅绘制数据有效范围内的等值线，自动忽略超出区域。

解决方法 #

方法1：标记无效值+指定原点对齐 #

1. 读取图像数据，将无效像素标记为NaN（根据数据实际情况调整过滤规则）：

import numpy as np

data = fits.open(path_sofia)[1].data[22]
# 示例：将低于3倍标准差的像素设为NaN，过滤背景噪声
data[data < 3 * data.std()] = np.nan

2. 调用imshow时指定origin='lower'，让像素坐标与WCS原点对齐，同时NaN值会被自动忽略：

a.imshow(data, transform=a.get_transform(wcs_sofia), origin='lower', cmap='viridis')

方法2：限制坐标轴范围为图像实际覆盖区域 #

先计算原图像在银道坐标系下的边界，再限制绘图范围：

# 获取原图像四个角的像素坐标
x_pix = [0, data.shape[1], data.shape[1], 0]
y_pix = [0, 0, data.shape[0], data.shape[0]]
# 转换为RA/DEC世界坐标
lon, lat = wcs_sofia.all_pix2world(x_pix, y_pix, 0)
# 再转换为银道坐标系的像素坐标
gal_lon, gal_lat = wcs_galactic.all_world2pix(lon, lat, 0)

# 设置坐标轴范围为图像实际覆盖的银道区域
a.set_xlim(min(gal_lon), max(gal_lon))
a.set_ylim(min(gal_lat), max(gal_lat))

# 绘制图像
a.imshow(data, transform=a.get_transform(wcs_sofia), origin='lower')

这样imshow仅绘制原图像在银道系内的有效区域，不会填充超出范围的冗余像素。

方法3：使用pcolormesh替代imshow #

pcolormesh基于网格顶点绘制，适配非笛卡尔坐标系的能力更强，不会产生冗余像素：

# 创建像素网格坐标
y, x = np.mgrid[0:data.shape[0], 0:data.shape[1]]
# 将像素坐标转换为银道世界坐标
gal_x, gal_y = wcs_sofia.all_pix2world(x, y, 0)
# 绘制图像
a.pcolormesh(gal_x, gal_y, data, transform=a.get_transform(wcs_galactic), cmap='viridis')

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