如何使用healpy基于HEALPix像素化绘制天区坐标2D直方图

我手里有一组天区天体的坐标数据，用下面的代码可以生成示例数据：

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

我现在想借助healpy工具，通过HEALPix像素化的方式，在Mollweide投影上绘制这些(l, b)坐标对应的点密度2D直方图，但折腾了好久都没搞定。

我查了healpy的官方教程，但教程里只讲了怎么绘制一维数组或者读取FITS文件绘图，完全没提到怎么用HEALPix像素化来实现这种点密度的2D直方图。

我自己尝试过调用hp.projaxes.MollweideAxes.hist2d(l, b, bins=10)这个函数，但根本运行不起来，现在彻底卡壳了。

不过目前我已经能用matplotlib在Mollweide投影上画出这些点了，代码是这样的：

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

有没有大佬能教教我怎么用healpy完成点密度2D直方图的绘制呀？

内容的提问来源于stack exchange，提问作者AlbertBranson