plot HEALPix maps in 3D using matplotlib

plot_hp3d.py

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
import matplotlib as mpl


def plot_hp3d(ax, m, r=1.0, vmin=None, vmax=None, cmap=None, norm=None,
              alpha=None, shade=False, **kwargs):
    m = np.asanyarray(m)
    *nmaps, npix = m.shape
    nside = hp.npix2nside(npix)

    r = np.broadcast_to(r, nmaps)

    if norm is None:
        norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
    cm = plt.cm.get_cmap(cmap)

    if alpha is not None:
        alpha = np.broadcast_to(alpha, nmaps)
    else:
        alpha = np.full(nmaps, None)

    out = np.empty(nmaps, dtype=object)
    for i in np.ndindex(*nmaps):
        pix = np.where(m[i] != hp.UNSEEN)[0]

        ver = np.empty((4*len(pix), 3))
        tri = np.empty((2*len(pix), 3), dtype=int)
        for j, p in enumerate(pix):
            ver[4*j:4*(j+1)] = r[i]*hp.boundaries(nside, p, 1).T
            tri[2*j:2*(j+1)] = [[4*j, 4*j+1, 4*j+3], [4*j+1, 4*j+2, 4*j+3]]

        sur = ax.plot_trisurf(ver[:, 0], ver[:, 1], ver[:, 2], triangles=tri,
                              alpha=alpha[i], shade=shade, antialiased=False,
                              **kwargs)

        col = cm(norm(m[i][pix]))

        sur.set_fc(np.repeat(col, 2, axis=0))
        sur.set_ec('none')
        sur.set_lw(0)

        out[i] = sur

    if not nmaps:
        out = out.item()

    return out


m = hp.read_map('cmb.fits')
m = hp.ud_grade(m, 64)

hp.mollview(m, min=-1e-4, max=1e-4, cmap='coolwarm')
plt.show()

ax = plt.subplot(111, projection='3d')
plot_hp3d(ax, m, vmin=-1e-4, vmax=1e-4, cmap='coolwarm', shade=True)
plt.show()

z_download.sh

curl -o cmb.fits 'http://pla.esac.esa.int/pla-sl/data-action?MAP.MAP_OID=13491'