In [10]:
from PIL import Image
import numpy as np
from matplotlib.image import pil_to_array
In [47]:
# Earth image extracted from basemap:
grayscale_pil_image ="./shadedrelief.jpg").convert("L")
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")
<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")
<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")
<matplotlib.text.Text at 0x292ef190>