def drawHealpixMap(hpxmap, lon, lat, size=1.0, xsize=501, coord='GC', **kwargs):
    """
    Draw local projection of healpix map.
    """
    ax = plt.gca()
    x = np.linspace(-size,size,xsize)
    y = np.linspace(-size,size,xsize)
    xx, yy = np.meshgrid(x,y)
    
    coord = coord.upper()
        
    if coord == 'GC':
        #Assumes map and (lon,lat) are Galactic, but plotting celestial
        llon, llat = image2sphere(*gal2cel(lon,lat),x=xx.flat,y=yy.flat)
        pix = ang2pix(get_nside(hpxmap),*cel2gal(llon,llat))
    elif coord == 'CG':
        #Assumes map and (lon,lat) are celestial, but plotting Galactic
        llon, llat = image2sphere(*cel2gal(lon,lat),x=xx.flat,y=yy.flat)
        pix = ang2pix(get_nside(hpxmap),*gal2cel(llon,llat))
    else:
        #Assumes plotting the native coordinates
        llon, llat = image2sphere(lon,lat,xx.flat,yy.flat)
        pix = ang2pix(get_nside(hpxmap),llon,llat)

    values = hpxmap[pix].reshape(xx.shape)
    zz = np.ma.array(values,mask=(values==hp.UNSEEN),fill_value=np.nan)

    return drawProjImage(xx,yy,zz,coord=coord,**kwargs)