Using Healpy to make star chart

Question

I'm using healpy to plot the locations of galaxies on the sky from a list of RAs and Decs. So far, I think I've been able to correctly plot the galaxies, but I'd like to improve the finished product. Is there any way to bin the number of galaxies that appear in each healpy tile, rather than just coloring being based on whether there is or isn't a catalog member in the tile?

Here I show the image that I'm currently making —

right now it's only really useful for telling you where the Milky Way isn't. Here's the code I'm using.

phis = [np.deg2rad(ra) for ra in ra_list]
thetas = [np.pi / 2 - np.deg2rad(dec) for dec in dec_list]

pixel_indices = hp.ang2pix(NSIDE, thetas, phis)

m = np.zeros(hp.nside2npix(NSIDE))
m[pixel_indices] = np.ones(num_galaxies_to_plot)


hp.mollview(m, title = 'Sky Locations of GLADE Galaxies', cbar = False, rot=(180, 0, 180), cmap = 'binary')
hp.graticule()

Answer 1

You could use numpy.bincount to create an array of the number of galaxies per pixels and then create a map of that.