How to create a HEALpix map with Python?

Question

I have a list of Galactic coordinates and want to plot an all-aky HEALpix map with it. Is there a library for that?

Answer 1

healpy has a function for that. Here is an example:

import healpy as hp
from healpy.newvisufunc import projview, newprojplot

import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np

def coords_to_hpx(lon, lat, nside, radec=False):
    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

    theta = np.radians(90. - b)
    phi = np.radians(l)
    indices = hp.ang2pix(nside, theta, phi)
    idx, counts = np.unique(indices, return_counts=True)
    hpx_map = np.zeros(npix, dtype=int)
    hpx_map[idx] = counts
    return hpx_map

data = np.loadtxt('photon_lists/p2209fb.txt')
l = data[:,3]
b = data[:,4]
l[l >= 180] -= 360

hpx = coords_to_hpx(l, b, nside=8, radec=False)
hpx[hpx > 50] = 50
projview(
        hpx,
        coord=["G"],
        graticule=True,
        graticule_labels=True,
        unit="counts",
        xlabel="longitude",
        ylabel="latitude",
        cb_orientation="horizontal",
        cbar_ticks=np.round(np.linspace(min(hpx), max(hpx), 10, endpoint=True)),
        projection_type="aitoff",
        )

plt.show()