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

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)

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")

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")

hp.gnomview(healpix_map, rot=(12.5, 41.9), reso=.5, xsize=1600, cmap="gray", flip="geo")
title("Gnomonic projection of Italy")