Make A 2d Histogram With Healpix Pixellization Using Healpy

The data are coordinates of objects in the sky, for example as follows: import pylab as plt import numpy as np l = np.random.uniform(-180, 180, 2000) b = np.random.uniform(-90, 90,

Solution 1:

Great question! I've written a short function to convert a catalogue into a HEALPix map of number counts:

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

defcat2hpx(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

You can then use that to populate the HEALPix map:

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

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

Here, the nside determines how fine or coarse your pixel grid is.

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

Also note that by sampling uniformly in Galactic latitude, you'll prefer data points at the Galactic poles. If you want to avoid that, you can scale that down with a cosine.

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