使用HEALPy节省像素数

Question

我有一个星系列表，可以将其绘制到healpy地图上(我正在使用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。

有没有一种方法可以在每个坐标占用多少像素之前计算出来，然后在此基础上改变倾倒的流量？

提前感谢！

Answer 1

函数hp.mollview的xsize参数只能用于绘图目的。如果要操作地图的分辨率，请使用hp.pixelfunc.ud_grade

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

在这行之后

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

在power=-2中使用ud_grade函数，因此总通量可以守恒：

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。