Healpy中正射投影像素值数量不符问题及像素面积咨询

问题描述 #

使用Healpy处理天图时遇到两个问题：

1. 在Mollweide投影（mollview）中为N个Healpix像素赋值100，但转换为正射投影（orthview）后，平面图像中值为100的像素数量与N不一致。
2. 询问Healpy中正射投影的像素面积是多少。

相关代码：

import healpy as hp
import numpy as np
import matplotlib.pyplot as plt

nside =64
npix = hp.nside2npix(nside)
hpxmap = np.arange(npix)

theta = np.pi/2. - np.deg2rad(0)
phi = np.deg2rad(35)
vec = hp.ang2vec(theta, phi, lonlat=False)
disc = hp.query_disc(nside, vec, radius=np.radians(10), nest=False)        
hpxmap[disc] = 100
hp.mollview(hpxmap, coord = ['C'], nest = False, cbar = True)
print(len(disc))


orth_proj = hp.orthview(hpxmap, half_sky=True,return_projected_map=True, rot=(0, 0), xsize= 100, norm = 'log', coord=[ 'C'])
orth_proj[np.where(orth_proj == orth_proj[0,0])] = 0
len(np.where(orth_proj == 100)[0][:])

像素数量不一致的原因与解决思路 #

核心原因 #

Healpix的球面像素是等面积划分，每个像素的球面面积固定（公式：4π / npix，其中npix=12*nside²）；而正射投影（orthview）是将球面投影到平面的非等面积透视投影，平面像素与球面像素不存在一对一的映射关系：

• 单个球面像素可能被投影到多个平面像素（尤其是靠近投影中心的区域）
• 多个边缘球面像素可能被合并到一个平面像素（投影边缘区域）
• 部分球面像素的边缘可能超出正射投影的可视范围，导致只被部分采样，不会被标记为100

另外，代码中设置的xsize=100决定了正射投影的平面分辨率，分辨率越低，采样误差越大，数量差异越明显。

解决思路 #

1. 避免直接对比平面与球面像素数：两者属于不同采样体系，没有严格的对应关系，不应以此作为验证标准。
2. 反向映射统计：如果需要统计原disc区域在正射投影中的覆盖，可遍历正射投影的每个平面像素，计算其对应的球面坐标，再查找对应的Healpix像素是否属于disc，以此统计数量。示例代码片段：

   # 获取正射投影的坐标网格
   theta, phi = hp.orthview(hpxmap, half_sky=True, return_projected_map=False, return_angles=True, rot=(0,0), xsize=100)
   # 将坐标转换为Healpix像素索引
   pix_indices = hp.ang2pix(nside, theta, phi, lonlat=False)
   # 统计属于disc的像素数量
   count = np.sum(np.isin(pix_indices, disc))
   

3. 提高投影分辨率：增大xsize参数（如设置为500或1000），可以减少采样带来的误差，让统计结果更接近真实的覆盖范围，但无法完全消除非等面积投影的固有差异。

正射投影像素面积说明 #

Healpy正射投影的平面像素没有固定的球面面积，因为正射投影的透视特性：

• 靠近投影中心（天顶）的平面像素，对应的球面区域面积较小
• 靠近投影边缘（地平线）的平面像素，对应的球面区域面积较大

如果需要计算某个平面像素对应的球面面积，可按以下步骤操作：

1. 获取该平面像素对应的球面坐标（theta, phi）
2. 计算该坐标附近微小区域的球面面积，或通过hp.pixelfunc.pixel_area获取对应Healpix像素的面积，再结合投影的重叠关系估算。

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