Representations of the sphere and the cosmic microwave background (CMB)¶
Autor: Eduardo Martín Calleja
In this post, which can be considered a continuation of the previous post, dedicated to the Mollweide projection, we will see an introduction to the Python healpy package. Healpy provides access from Python to the HEALPix set of functions, which are the standard for data representation from the various missions that measure the temperature of the microwave background radiation of the universe.
HALPIx provides an algorithm for subdividing the surface of the sphere in a series of picture elements (pixels) with the property that they all represent exactly the same area of the original spherical surface, and furthermore, these pixels are arranged into lines of equal latitude. Then this subdivision can be transferred to a plane projection like Mollweide, as it preserve the equal-area property in the flat representation.
Having graphic representations of the spherical surface with pixels of the same area is crucial for the presentation and subsequent analysis of the distribution of mass, energy, radiation, etc.. Allowing relative densities to be compared and to apply various algorithms. In particular, the library HEALPix was originally developed to represent the distribution of the microwave background radiation which has been compiling data on various missions, although its use has since spread to other fields.
Imports and references¶
%matplotlib inline
from __future__ import division
import numpy as np
import healpy as hp
import astroML
#This removes some nasty deprecation warnings that do not interfere with the execution
import warnings
warnings.filterwarnings('ignore')
# This IPython magic generates a table with version information
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information numpy, healpy, astroML, astroPy
Some references:
Some additional information can be found in the HEALPix site
Coordinates and pixel numbering¶
In healpy spherical coordinates are used to position the pixels: the polar angle "theta" $ $ measured from the North pole (colatitude), which takes values in the interval $ [0, ] \(, and the azimuth angle "phi" \) $ (longitude) that takes values in the range $ [0,2 ] $, measured to the East of the globe.
The resolution of the subdivision in pixels is defined by parameter NSIDE
The function * nside2npix (NSIDE) * returns the number of pixels that will decompose the spherical surface for a given NSIDE. The minimum resolution is 12 pixels, and every time you double NSIDE (which should always be a power of 2) each pixel is divided into four, as shown in the following table:
for NSIDE in 2.**np.arange(6):
print 'The number of pixels for NSIDE = %2d is %5d' %(NSIDE, hp.nside2npix(NSIDE))
The pixels are numbered by default according to the scheme RING, consisting in enumerating them consecutively along rings of constant latitude, starting by the northern-most ring of pixels.
The pix2ang(nside, ipix) function returns the angular coordinates for a given pixel number of array of pixel numbers. Note that \(\theta\) takes values between 0 and \(\pi\) radians, and \(\phi\) between 0 and $ 2 $ radians. Let's see an example:
NSIDE = 1
ipix = np.arange(hp.nside2npix(NSIDE))
print 'The numbering of the pixels for NSIDE = %d is: \n' %NSIDE, ipix
theta, phi = hp.pix2ang(NSIDE,ipix)
print '\n array theta = ', theta/np.pi, 'radians'
print '\n array phi = ', phi/np.pi, 'radians'
Conversely, the function ang2pix(nside, theta, phi) returns the numbering of pixels for a given set of angular coordinates (theta, phi), which as before can be single or two separate coordinate arrays
NSIDE = 1
hp.ang2pix(NSIDE, np.pi/2, np.pi)
Let's see what happens now if we use two separate arrays with the angular coordinates:
# If arrays of the same length: the matching is element to element
theta = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4])
phi = np.array([0, np.pi/2, np.pi, 3*np.pi/2])
hp.ang2pix(NSIDE, theta, phi)
# Case of arrays of different lengths: we get an exception
theta = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4])
phi2 = np.array([0, np.pi/2, np.pi, 3*np.pi/2, 3.5*np.pi/2])
try:
hp.ang2pix(NSIDE, theta, phi2)
except Exception as e:
print "Atttention: ",e
# But if we change the theta array to a column array
# Then the numpy broadcasting is possible,
# and returns a 2D array of pixel numbers
theta2 = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4]).reshape(-1,1)
phi2 = np.array([0, np.pi/2, np.pi, 3*np.pi/2, 3.5*np.pi/2])
hp.ang2pix(NSIDE, theta2, phi2)
Graphical representation in Mollweide projection¶
The healpy function to plot an array of pixels in the Mollweide projection is mollview(). To be represented with mollview() the array must be one-dimensional, and its length must match nside2npix(NSIDE). Each array element will then correspond to a pixel in the HEALPix subdivision. Here's an example:
npixels = hp.nside2npix(32)
np.random.seed(0)
test = np.random.rand(npixels)
hp.mollview(test, min=0, max=1,
title = 'Mollweide projection test',
unit='This is a description of the measurement unit')
We can see the pixels numbering scheme in the following example, where at the same time we will introduce the functions projscatter() and projtext(). It will use the default schema: "RING":
NSIDE = 2
npixels = hp.nside2npix(NSIDE)
img = np.linspace(0, 255, num=npixels)
index = np.arange(npixels)
theta, phi= hp.pix2ang(NSIDE,index)
hp.mollview(img, min=0, max = 255, unit='Range of values')
hp.projscatter(theta, phi)
for i in index:
hp.projtext(theta[i]-0.05, phi[i], i)
Note that in the above scheme, the rendering order of the pixels is above (North) to below (South) and right (West) to left (East), as corresponds to a projection of the sky in which the East is located to the left. If you wish to represent a geographical map, with the East to the right, you must use the flip = "geo" option.
healpy works internally with angle coordinates theta and phi in radians, but you can also provide coordinates in other systems, such as galactic coordinates (lon, lat) or equatorial coordinates (RA, Dec). For this reason, the functions projscatter() and projtext() support the option lonlat = True, indicating that the coordinates we provide are not theta (colatitude) and phi (longitude) in radians, but a longitude and latitude in degrees. For example, in equatorial coordinates the angular coordinates would be right ascension and declination.
Let's see an example following this post
Again, recall that the East is to the left of the figure
lon = [30, 105, 270]
lat = [45, -30, 60]
hp.mollview(title="Mollweide in (lon, lat) coordinates")
hp.graticule() # add grid
hp.projscatter(lon, lat, lonlat=True)
hp.projtext(30, 45,'(30,45)', lonlat=True)
hp.projtext(105, -30,'(105,-30)', lonlat=True)
hp.projtext(270, 60,'(270,60)', lonlat=True);
Download and display of data from the Cosmic Microwave Background (CMB) WMAP mission¶
To access the WMAP data of the cosmic microwave background we will use the astroML Python module. The function to get an array ready for graphing with healpy is:
astroML.datasets.fetch_wmap_temperatures()
You can download two versions of this data: A "full" version that includes the pollutant frequencies, mainly radiation from our own galaxy or other extragalactic sources, or a version in which the frequency bands which concentrates this contamination have been masked.
Let's download the unmasked data:
from astroML.datasets import fetch_wmap_temperatures
wmap_unmasked = fetch_wmap_temperatures(masked=False)
wmap_unmasked.shape
These data are in the form of an array of Numpy already prepared for plotting with healpy. It can be seen that the level of graphic resolution is that corresponding to a NSIDE of 512
hp.nside2npix(512)
Let us now proceed to its graphical representation:
hp.mollview(wmap_unmasked, min=-1, max=1,
title='Unmasked WMAP data \n in galactic coordinates',
unit=r'$\Delta$T (mK)', xsize = 200)
The ordering of the pixels in the array above corresponds to the galactic coordinate system. We can see the presence in the chart, at the position lat = 0, of a band of thermal pollution from our galaxy. This is irrelevant when it comes to construct the healpy graphical representation, since all it does is reading the array and place each pixel in the position that corresponds to the numbering scheme adopted (RING default), and to the resolution indicated by the length of the array.
However, we can tell healpy that we are working on a specific coordinate system and want to change it to another different coordinate system. healpy can work with three types of coordinates: "G": Galactic, "E": Ecliptic and "C": Celestial Equator. We can for example redraw the graph, converting galactic to equatorial coordinates by using the option: coord = 'GC'.
Furthermore, we can even rotate the chart above to center the origin in, for example, (RA, dec) = (180.0), thus obtaining a map that shows the status of the "Milky Way" in equatorial coordinates centered at RA = 180 degrees
hp.mollview(wmap_unmasked, min=-1, max=1,
title='Unmasked WMAP data \n in equatorial coordinates',
coord='GC', rot=(180,0),
unit=r'$\Delta$T (mK)', xsize = 200)
hp.graticule()
hp.projtext(180, 2,'RA=180', lonlat=True)
hp.projtext(90, 2, 'RA=90', lonlat=True)
hp.projtext(270, 2, 'RA=270', lonlat=True);
No hay comentarios:
Publicar un comentario