使用Astropy从HDF5文件绘制法拉第天空均值与标准差天图时的纬度范围错误修复求助

你的错误核心在于：把法拉德深度的均值/标准差数值当成了银河坐标（l/b）来使用，这完全搞反了数据含义！均值/std的数值是天区像素的测量结果，不是坐标值，这才导致了纬度超出-90°~90°范围的报错。

下面是具体的修复步骤和完整代码：

1. 先理清数据结构 #

Faraday Sky 2020的HDF5文件中，faraday_sky_mean和faraday_sky_std是2D像素数组（维度为(360, 720)）：

• 第一维度对应银河纬度b，范围是-90°到90°，共360个采样点
• 第二维度对应银河经度l，范围是0°到360°，共720个采样点
  数组内的数值是对应天区的法拉德深度均值/标准差，而非坐标。

2. 生成正确的银河坐标网格 #

我们需要先创建与像素数据匹配的银河坐标网格：

# 根据数据维度生成坐标序列
l = np.linspace(0, 360, faraday_sky_mean.shape[1]) * u.degree
b = np.linspace(-90, 90, faraday_sky_mean.shape[0]) * u.degree
# 生成网格坐标（每个像素对应一对l/b）
l_grid, b_grid = np.meshgrid(l, b)

3. 适配Aitoff投影的坐标转换 #

Aitoff投影要求经度范围是-180°到180°，所以需要把经度l wrap到这个区间后转成弧度：

c = SkyCoord(l=l_grid, b=b_grid, frame='galactic')
l_rad = c.l.wrap_at(180*u.deg).radian
b_rad = c.b.radian

4. 用2D绘图展示数据 #

之前的plt.plot是画散点，对于2D天图，应该用plt.pcolormesh通过颜色映射展示数据分布：

完整修复后的代码 #

from astropy import units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
import h5py

with h5py.File('faraday2020.hdf5','r') as hdf:
    print(list(hdf.keys()))
    faraday_sky_mean = hdf['faraday_sky_mean'][:]
    faraday_sky_std = hdf['faraday_sky_std'][:]
    # 打印数据基本信息，确认数据范围
    print(f"Mean shape: {faraday_sky_mean.shape}, dtype: {faraday_sky_mean.dtype}")
    print(f"Mean range: {faraday_sky_mean.min()} ~ {faraday_sky_mean.max()} rad/m²")
    print(f"Std shape: {faraday_sky_std.shape}, dtype: {faraday_sky_std.dtype}")
    print(f"Std range: {faraday_sky_std.min()} ~ {faraday_sky_std.max()} rad/m²")

# 生成匹配数据的银河坐标网格
l = np.linspace(0, 360, faraday_sky_mean.shape[1]) * u.degree
b = np.linspace(-90, 90, faraday_sky_mean.shape[0]) * u.degree
l_grid, b_grid = np.meshgrid(l, b)

# 转换为适配Aitoff投影的弧度坐标
c = SkyCoord(l=l_grid, b=b_grid, frame='galactic')
l_rad = c.l.wrap_at(180*u.deg).radian
b_rad = c.b.radian

# 创建子图分别展示均值和标准差
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6), subplot_kw={'projection': 'aitoff'})

# 绘制均值天图
ax1.set_title("Faraday Depth Mean", y=1.08, fontsize=12)
ax1.grid(True)
im1 = ax1.pcolormesh(l_rad, b_rad, faraday_sky_mean, cmap='viridis', shading='auto')
fig.colorbar(im1, ax=ax1, label='Mean Faraday Depth (rad/m²)')

# 绘制标准差天图
ax2.set_title("Faraday Depth Standard Deviation", y=1.08, fontsize=12)
ax2.grid(True)
im2 = ax2.pcolormesh(l_rad, b_rad, faraday_sky_std, cmap='plasma', shading='auto')
fig.colorbar(im2, ax=ax2, label='Std Dev Faraday Depth (rad/m²)')

# 调整布局避免标题重叠
plt.subplots_adjust(top=0.9, bottom=0.1, wspace=0.3)
plt.show()

关键修复点总结 #

• 停止将均值/std数值错误作为坐标，改为生成匹配数据的银河坐标网格
• 使用pcolormesh绘制2D天图，通过颜色映射直观展示天区数据分布
• 正确处理Aitoff投影的经度范围（wrap到-180°~180°）

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