3d_map_planck_matplotlib.py

#from mayavi import mlab
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt

# load a Planck map, Kcmb -> mKcmb
m = hp.read_map("HFI_SkyMap_100_2048_R3.01_full.fits", 0) * 1e3
nside = hp.npix2nside(len(m))

vmin = -1; vmax = 1

# size of the grid
xsize = ysize = 500

theta = np.linspace(np.pi, 0, ysize)
phi   = np.linspace(-np.pi, np.pi, xsize)
longitude = np.radians(np.linspace(-180, 180, xsize))
latitude = np.radians(np.linspace(-90, 90, ysize))

# project the map to a rectangular matrix xsize x ysize
PHI, THETA = np.meshgrid(phi, theta)
grid_pix = hp.ang2pix(nside, THETA, PHI)
grid_map = m[grid_pix]

# Create a sphere
r = 1
x = r*np.sin(THETA)*np.cos(PHI)
y = r*np.sin(THETA)*np.sin(PHI)
z = r*np.cos(THETA)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
from matplotlib import cm

grid_map = np.clip(grid_map, -1, 1)
grid_map -= grid_map.min()
grid_map /= grid_map.max()

ax.plot_surface( x, y,  z, cstride=1, rstride=1, facecolors=cm.jet(grid_map) )

# Set an equal aspect ratio
ax.set_aspect('equal')

plt.show()