使用 HEALPy 保存像素数答案

【问题标题】：Conserving pixel counts with HEALPy使用 HEALPy 保存像素数
【发布时间】：2017-07-26 14:13:14
【问题描述】：

我有一个星系列表要绘制到 healpix 地图上（我正在使用 healpy 来做）每个星系都有一个固定的通量，我需要让它们以这样的方式绘制，即每个星系的通量都是守恒的在地图上。

这是我的代码：

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

nside = 8
xsize = 100

ra = np.array([pi/4,pi/3])
dec = np.array([pi/4,pi/3])
flux = np.array([10,20])

hpm = np.zeros(hp.nside2npix(nside)) #Blank healpix map
pixindex = hp.ang2pix(nside, dec, ra)

np.add.at(hpm,pixindex,flux) #Add flux onto correct pixels

img=hp.mollview(hpm,coord=['E'],xsize=xsize,return_projected_map=True)
print(np.sum(img[img>0]))

我得到的结果是 140 而不是 30，这是通量的真正总和。

我知道发生了什么，并且相同的通量分布在多个像素上（第一个星系为 6 个像素，第二个星系为 4 个像素），我知道我可以这样做：

newimg = img * (np.sum(flux)/np.sum(img[img>0]))

这会保存总光子数，但不一定会保存我需要的每个星系的光子数。即这种方法最终导致第一个星系贡献 12.86 的通量，第二个星系贡献 17.14 的通量。

有没有办法在每个坐标占用多少像素之前计算出，然后根据此改变倾倒的通量量？

提前致谢！

【问题讨论】：

标签： python maps astronomy healpy

【解决方案1】：

函数hp.mollview 的xsize 参数应仅用于绘图目的。如果您想控制地图的分辨率，请使用hp.pixelfunc.ud_grade

例如，如果你想从nside=8 转到nside=32，

行后

np.add.at(hpm, pixindex, flux) #Add flux onto correct pixels

将ud_grade 函数与power=-2 一起使用，这样总通量可以守恒：

hpm_nside_32 = hp.pixelfunc.ud_grade(hpm, power=-2, nside_out=32)

总和

np.sum(hpm_nside_32)

将保存在30。

如果您需要 mollview 来保存助焊剂，我无法提供解决方案。我能得到的最接近的方法是根据 mollview 图像中的像素数与 hpm 中的像素数之比缩小 img 值。第一项 xsize * xsize / 2. 是 mollview 中的像素数，第二项 2. * np.pi / 8. 是具有半短轴长度一半长半轴 pi * r * 2r 的椭圆的面积与面积之比一个矩形4r * 2r。

len(hpm) / ((xsize * xsize / 2.) * (2. * np.pi / 8.)) * np.sum(img[img > 0])

当xsize = 100时，总和将变为27.380；当xsize = 1000 时，近似值要好得多，在29.981；当xsize = 10000 时，总通量变为29.988。

另一种为您提供良好近似值的方法是计算img 中非-inf 像素数与地图中像素数之比（nside=8 为768）：

float(len(hpm)) / float(np.sum(img>=0)) * np.sum(img[img>0])

在xsize = 100, 1000, 10000，通量将分别在28.295, 30.075, 29.998。

【讨论】：