使用healpy进行HEALPix像素化的2D直方图制作

Question

使用healpy进行HEALPix像素化的2D直方图制作

10

这些数据是天空中物体的坐标，例如如下：

import pylab as plt
import numpy as np
l = np.random.uniform(-180, 180, 2000)
b = np.random.uniform(-90, 90, 2000)

我想做一个二维直方图，以便在天空中使用HEALPix像素化和Mollweide投影绘制一些点的密度地图，这些点具有(l, b)坐标。我该如何使用healpy实现？

教程：

http://healpy.readthedocs.io/en/v1.9.0/tutorial.html

说明如何绘制1D数组或fits文件，但我找不到如何使用该像素化绘制2D直方图。

我也找到了这个函数，但它没有起作用，所以我被卡住了。

hp.projaxes.MollweideAxes.hist2d(l, b, bins=10)

我可以使用莫尔韦德投影将这些点绘制成图表：

l_axis_name ='Latitude l (deg)'
b_axis_name = 'Longitude b (deg)'

fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(111, projection="mollweide")
ax.grid(True)

ax.scatter(np.array(l)*np.pi/180., np.array(b)*np.pi/180.)

plt.show()

非常感谢您的帮助。

- AlbertBranson

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Daniel Lenz · Accepted Answer

非常好的问题！我写了一个简短的函数，将目录转换成HEALPix数量计数地图：

from astropy.coordinates import SkyCoord
import healpy as hp
import numpy as np

def cat2hpx(lon, lat, nside, radec=True):
    """
    Convert a catalogue to a HEALPix map of number counts per resolution
    element.

    Parameters
    ----------
    lon, lat : (ndarray, ndarray)
        Coordinates of the sources in degree. If radec=True, assume input is in the icrs
        coordinate system. Otherwise assume input is glon, glat

    nside : int
        HEALPix nside of the target map

    radec : bool
        Switch between R.A./Dec and glon/glat as input coordinate system.

    Return
    ------
    hpx_map : ndarray
        HEALPix map of the catalogue number counts in Galactic coordinates

    """

    npix = hp.nside2npix(nside)

    if radec:
        eq = SkyCoord(lon, lat, 'icrs', unit='deg')
        l, b = eq.galactic.l.value, eq.galactic.b.value
    else:
        l, b = lon, lat

    # conver to theta, phi
    theta = np.radians(90. - b)
    phi = np.radians(l)

    # convert to HEALPix indices
    indices = hp.ang2pix(nside, theta, phi)

    idx, counts = np.unique(indices, return_counts=True)

    # fill the fullsky map
    hpx_map = np.zeros(npix, dtype=int)
    hpx_map[idx] = counts

    return hpx_map

然后您可以使用此内容填充HEALPix地图：

l = np.random.uniform(-180, 180, 20000)
b = np.random.uniform(-90, 90, 20000)

hpx_map = hpx.cat2hpx(l, b, nside=32, radec=False)

在这里，nside 决定了您的像素网格是多么精细或粗糙。

hp.mollview(np.log10(hpx_map+1))

还要注意的是，如果按银道纬度均匀采样，您会更倾向于选择银河极点的数据点。如果想避免这种情况，可以通过余弦函数进行缩放。

hp.orthview(np.log10(hpx_map+1), rot=[0, 90])
hp.graticule(color='white')