from PIL import Image
import numpy as np
from matplotlib.image import pil_to_array
# Earth image extracted from basemap:
# https://github.com/matplotlib/basemap/blob/master/lib/mpl_toolkits/basemap/data/shadedrelief.jpg
grayscale_pil_image = Image.open("./shadedrelief.jpg").convert("L")
image_array = pil_to_array(grayscale_pil_image)
print image_array.shape
(5400, 10800)
import healpy as hp
theta = np.linspace(0, np.pi, num=image_array.shape[0])[:, None]
phi = np.linspace(-np.pi, np.pi, num=image_array.shape[1])
nside = 512
print "Pixel area: %.2f square degrees" % hp.nside2pixarea(nside, degrees=True)
Pixel area: 0.01 square degrees
pix = hp.ang2pix(nside, theta, phi)
healpix_map = np.zeros(hp.nside2npix(nside), dtype=np.double)
healpix_map[pix] = image_array
hp.mollview(healpix_map, cmap="gray", xsize=2000, flip="geo")
title("Mollweide view of the Earth")
<matplotlib.text.Text at 0x207bf090>
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")
<matplotlib.text.Text at 0x1f51a3d0>
hp.gnomview(healpix_map, rot=(12.5, 41.9), reso=.5, xsize=1600, cmap="gray", flip="geo")
title("Gnomonic projection of Italy")
<matplotlib.text.Text at 0x292ef190>