In [10]:
from PIL import Image
import numpy as np
from matplotlib.image import pil_to_array

In [47]:
# Earth image extracted from basemap:
image_array = pil_to_array(grayscale_pil_image)

In [48]:
print image_array.shape

(5400, 10800)

In [49]:
import healpy as hp

In [83]:
theta = np.linspace(0, np.pi, num=image_array.shape[0])[:, None]
phi = np.linspace(-np.pi, np.pi, num=image_array.shape[1])

In [97]:
nside = 512
print "Pixel area: %.2f square degrees" % hp.nside2pixarea(nside, degrees=True)

Pixel area: 0.01 square degrees

In [98]:
pix = hp.ang2pix(nside, theta, phi)

In [99]:
healpix_map = np.zeros(hp.nside2npix(nside), dtype=np.double)

In [100]:
healpix_map[pix] = image_array

In [101]:
hp.mollview(healpix_map, cmap="gray", xsize=2000, flip="geo")
title("Mollweide view of the Earth")

Out[101]:
<matplotlib.text.Text at 0x207bf090>
In [104]:
hp.mollview(hp.smoothing(healpix_map, fwhm=np.radians(2)), cmap="gray", xsize=2000, flip="geo")
title("Earth map smoothed with a 1 degree gaussian beam")

Out[104]:
<matplotlib.text.Text at 0x1f51a3d0>
In [109]:
hp.gnomview(healpix_map, rot=(12.5, 41.9), reso=.5, xsize=1600, cmap="gray", flip="geo")
title("Gnomonic projection of Italy")

Out[109]:
<matplotlib.text.Text at 0x292ef190>