基于Healpy的CMB地图南北纬10度像素条提取及角度问题

你完全不需要旋转坐标系！问题的核心是搞清楚Healpy球坐标系里的theta和地理纬度的对应关系，直接计算对应角度就能提取目标像素。

关键坐标系对应关系 #

Healpy使用的球坐标系参数定义：

• theta：极角，从**北极点（地理纬度+90°）**开始测量，范围是0到π弧度（对应0°到180°），终点是南极点（地理纬度-90°）
• phi：方位角，从0到2π弧度（对应0°到360°），就是常规的经度

地理纬度lat（北纬为正，南纬为负）和theta的转换公式非常简单：

theta = np.deg2rad(90 - lat)

举个例子：

• 北半球+10°纬度（赤道以北10°）：lat=10° → theta = 90° - 10° = 80° → 转弧度就是np.deg2rad(80)
• 南半球-10°纬度（赤道以南10°）：lat=-10° → theta = 90° - (-10°) = 100° → 转弧度就是np.deg2rad(100)

修正后的完整代码 #

下面是提取这两条纬度线像素的完整代码，同时纠正了你之前代码里的theta计算误区：

import numpy as np
import healpy as hp

# 读取CMB地图
fname = 'COM_CMB_IQU-070-fgsub-sevem-field-Pol_1024_R2.01_full.fits'
tmap = hp.read_map(fname)
nside = hp.get_nside(tmap)

# 定义目标纬度（北纬+10°，南纬-10°）
target_lat_north = 10
target_lat_south = -10

# 转换为Healpy的极角theta
theta_north = np.deg2rad(90 - target_lat_north)
theta_south = np.deg2rad(90 - target_lat_south)

# 生成方位角phi的范围（0到2π，覆盖整个经度）
phi_range = np.linspace(0, 2*np.pi, 1000)  # 这里的1000可以根据需要调整采样密度

# 提取北半球+10°纬度的像素索引
north_pix_indices = hp.ang2pix(nside, theta_north, phi_range)
# 提取南半球-10°纬度的像素索引
south_pix_indices = hp.ang2pix(nside, theta_south, phi_range)

# 可选：获取对应像素的CMB温度值
north_temperatures = tmap[north_pix_indices]
south_temperatures = tmap[south_pix_indices]

# 输出示例
print("北半球+10°纬度像素索引示例：", north_pix_indices[:5])
print("南半球-10°纬度像素索引示例：", south_pix_indices[:5])

补充说明 #

1. 为什么不用旋转？因为Healpy的坐标系已经直接覆盖了整个天球，只要正确计算目标纬度对应的theta，就能直接定位到南半球的纬线，完全不需要额外旋转坐标系。
2. 关于phi_range：用np.linspace生成足够多的方位角采样点，可以覆盖整条纬线上的所有像素（Healpy会自动匹配最近的像素索引）。采样点数量越多，越能完整覆盖这条纬线的像素。
3. 如果你之前的代码是想提取北极点以南10°的圈（也就是lat=80°），那你的theta计算是对的，但如果是赤道以北10°，就需要用上面的转换公式。

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