How to use the healpy.pixelfunc.ang2pix function in healpy

npix = (int(wcs.wcs.crpix[0] * 2), int(wcs.wcs.crpix[1] * 2))
    pix_crds = np.dstack(np.meshgrid(np.arange(npix[0]),
                                     np.arange(npix[1]))).swapaxes(0, 1).reshape((npix[0] * npix[1], 2))
    if wcs.wcs.naxis == 2:
        sky_crds = wcs.wcs_pix2world(pix_crds, 0)
    else:
        use_wcs = wcs.dropaxis(2)
        sky_crds = use_wcs.wcs_pix2world(pix_crds, 0)
        
    sky_crds *= np.radians(1.)
    sky_crds[0:, 1] = (np.pi / 2) - sky_crds[0:, 1]

    fullmask = np.isnan(sky_crds)
    mask = (fullmask[0:, 0] + fullmask[0:, 1]) == 0
    ipixs = -1 * np.ones(npix, int).T.flatten()
    ipixs[mask] = hp.pixelfunc.ang2pix(hpx.nside, sky_crds[0:, 1][mask],
                                       sky_crds[0:, 0][mask], hpx.nest)

    # Here we are counting the number of HEALPix pixels each WCS pixel points to;
    # this could probably be vectorized by filling a histogram.
    d_count = {}
    for ipix in ipixs:
        if ipix in d_count:
            d_count[ipix] += 1
        else:
            d_count[ipix] = 1

    # Here we are getting a multiplicative factor that tells use how to split up
    # the counts in each HEALPix pixel (by dividing the corresponding WCS pixels
    # by the number of associated HEALPix pixels).
    # This could also likely be vectorized.
    mult_val = np.ones(ipixs.shape)

nest (Optional[:obj:`bool`]): If :obj:`True` (the default), nested pixel ordering
            will be used. If :obj:`False`, ring ordering will be used.

    Returns:
        The HEALPix pixel index or indices. Has the same shape as the input :obj:`l`
        and :obj:`b`.
    """

    theta = np.radians(90. - b)
    phi = np.radians(l)

    if not hasattr(l, '__len__'):
        if (b < -90.) or (b > 90.):
            return -1

        pix_idx = hp.pixelfunc.ang2pix(nside, theta, phi, nest=nest)

        return pix_idx

    idx = (b >= -90.) & (b <= 90.)

    pix_idx = np.empty(l.shape, dtype='i8')
    pix_idx[idx] = hp.pixelfunc.ang2pix(nside, theta[idx], phi[idx], nest=nest)
    pix_idx[~idx] = -1

    return pix_idx

theta = np.radians(90. - b)
    phi = np.radians(l)

    if not hasattr(l, '__len__'):
        if (b < -90.) or (b > 90.):
            return -1

        pix_idx = hp.pixelfunc.ang2pix(nside, theta, phi, nest=nest)

        return pix_idx

    idx = (b >= -90.) & (b <= 90.)

    pix_idx = np.empty(l.shape, dtype='i8')
    pix_idx[idx] = hp.pixelfunc.ang2pix(nside, theta[idx], phi[idx], nest=nest)
    pix_idx[~idx] = -1

    return pix_idx

def _lbIndx(self,l,b):
        """Return the index in the _combineddata array corresponding to this (l,b)"""
        for nside in self._nsides:
            # Search for the pixel in this Nside level
            tpix= healpy.pixelfunc.ang2pix(nside,(90.-b)*_DEGTORAD,
                                           l*_DEGTORAD,nest=True)
            indx= (self._pix_info['healpix_index'] == tpix)\
                *(self._pix_info['nside'] == nside)
            if numpy.sum(indx) == 1:
                return self._indexArray[indx]
            elif numpy.sum(indx) > 1:
                raise IndexError("Given (l,b) pair has multiple matches!")
        raise IndexError("Given (l,b) pair not within the region covered by the extinction map")

# convert RA and Dec to phi and theta coordinates
phi = ra * np.pi/180.
theta = (90.0 - dec) * np.pi/180.

# to just plot all points, do
hp.visufunc.projplot(theta, phi, 'bo')
hp.visufunc.graticule()
# more examples are here
# https://healpy.readthedocs.org/en/latest/generated/healpy.visufunc.projplot.html#healpy.visufunc.projplot

## to plot a 2D histogram in the Mollweide projection
# define the HEALPIX level
NSIDE = 32

# find the pixel ID for each point
pix = hp.pixelfunc.ang2pix(NSIDE, theta, phi)

# select all points above Dec > -30
pix_unique = np.unique(pix[dec > -30])

# prepare the map array
m = np.zeros(hp.nside2npix(NSIDE), dtype='u2')

# tag the map array with pixels above Dec > -30
m[pix_unique] = 1

# plot map ('C' means the input coordinates were in the equatorial system)
hp.visufunc.mollview(m, coord=['C'], rot=(0, 0, 0), notext=True, title='', cbar=False)
hp.visufunc.graticule()

plt.show()

# to just plot all points, do
#hp.visufunc.projplot(theta, phi, 'bo')
#hp.visufunc.projplot(theta_lamost, phi_lamost, 'bo')
#hp.visufunc.graticule() # just the bare background w/ lines
# more examples are here
# https://healpy.readthedocs.org/en/latest/generated/healpy.visufunc.projplot.html#healpy.visufunc.projplot

## to plot a 2D histogram in the Mollweide projection
# define resolution of map
# NSIDE = 32 
NSIDE =  128

# find the pixel ID for each point
# pix = hp.pixelfunc.ang2pix(NSIDE, theta, phi)
pix_lamost = hp.pixelfunc.ang2pix(NSIDE, theta_lamost, phi_lamost)
pix_apogee = hp.pixelfunc.ang2pix(NSIDE, theta_apogee, phi_apogee)
pix_all = hp.pixelfunc.ang2pix(NSIDE, theta_all, phi_all)
# pix is in the order of ra and dec

# prepare the map array
m_lamost = hp.ma(np.zeros(hp.nside2npix(NSIDE), dtype='float'))
mask_lamost = np.zeros(hp.nside2npix(NSIDE), dtype='bool')

for pix_val in np.unique(pix_lamost):
    choose = np.where(pix_lamost==pix_val)[0]
    if len(choose) == 1:
        m_lamost[pix_val] = am_lamost[choose[0]]
    else:
        m_lamost[pix_val] = np.median(am_lamost[choose])

mask_lamost[np.setdiff1d(np.arange(len(m_lamost)), pix_lamost)] = 1
m_lamost.mask = mask_lamost

#hp.visufunc.projplot(theta_lamost, phi_lamost, 'bo')
#hp.visufunc.graticule() # just the bare background w/ lines
# more examples are here
# https://healpy.readthedocs.org/en/latest/generated/healpy.visufunc.projplot.html#healpy.visufunc.projplot

## to plot a 2D histogram in the Mollweide projection
# define the HEALPIX level
# NSIDE = 32 # defines the resolution of the map
# NSIDE =  128 # from paper 1
NSIDE = 64

# find the pixel ID for each point
# pix = hp.pixelfunc.ang2pix(NSIDE, theta, phi)
pix_lamost = hp.pixelfunc.ang2pix(NSIDE, theta_lamost, phi_lamost)
pix_apogee = hp.pixelfunc.ang2pix(NSIDE, theta_apogee, phi_apogee)
pix_all = hp.pixelfunc.ang2pix(NSIDE, theta_all, phi_all)
# pix is in the order of ra and dec

# prepare the map array
m_lamost = hp.ma(np.zeros(hp.nside2npix(NSIDE), dtype='float'))
mask_lamost = np.zeros(hp.nside2npix(NSIDE), dtype='bool')

for pix_val in np.unique(pix_lamost):
    choose = np.where(pix_lamost==pix_val)[0]
    if len(choose) == 1:
#         #m_lamost[pix_val] = rmag_lamost[choose[0]]
        m_lamost[pix_val] = val_lamost[choose[0]]
    else:
        #m_lamost[pix_val] = np.median(rmag_lamost[choose])
        m_lamost[pix_val] = np.median(val_lamost[choose])

mask_lamost[np.setdiff1d(np.arange(len(m_lamost)), pix_lamost)] = 1

# to just plot all points, do
#hp.visufunc.projplot(theta, phi, 'bo')
#hp.visufunc.projplot(theta_lamost, phi_lamost, 'bo')
#hp.visufunc.graticule() # just the bare background w/ lines
# more examples are here
# https://healpy.readthedocs.org/en/latest/generated/healpy.visufunc.projplot.html#healpy.visufunc.projplot

## to plot a 2D histogram in the Mollweide projection
# define the HEALPIX level
# NSIDE = 32 # defines the resolution of the map
# NSIDE =  128 # from paper 1
NSIDE = 64

# find the pixel ID for each point
# pix = hp.pixelfunc.ang2pix(NSIDE, theta, phi)
pix_lamost = hp.pixelfunc.ang2pix(NSIDE, theta_lamost, phi_lamost)
pix_apogee = hp.pixelfunc.ang2pix(NSIDE, theta_apogee, phi_apogee)
pix_all = hp.pixelfunc.ang2pix(NSIDE, theta_all, phi_all)
# pix is in the order of ra and dec

# prepare the map array
m_lamost = hp.ma(np.zeros(hp.nside2npix(NSIDE), dtype='float'))
mask_lamost = np.zeros(hp.nside2npix(NSIDE), dtype='bool')

for pix_val in np.unique(pix_lamost):
    choose = np.where(pix_lamost==pix_val)[0]
    if len(choose) == 1:
#         #m_lamost[pix_val] = rmag_lamost[choose[0]]
        m_lamost[pix_val] = val_lamost[choose[0]]
    else:
        #m_lamost[pix_val] = np.median(rmag_lamost[choose])
        m_lamost[pix_val] = np.median(val_lamost[choose])