1

我有一个星系列表要绘制到 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 的通量。

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

提前致谢！

4

1 回答 1

1

该xsize函数的参数hp.mollview应仅用于绘图目的。如果您想操纵地图的分辨率，请使用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。

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

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

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

于 2017-08-21T17:31:52.613 回答