-->

Etiquetas

Representations of the sphere and the cosmic microwave background (CMB)

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

In [1]:
%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
Out[1]:
SoftwareVersion
Python2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython2.3.1
OSLinux 3.13.0 45 generic x86_64 with debian jessie sid
numpy1.9.1
healpy1.8.4
astroML0.2
astroPy0.4.3
Fri Feb 20 21:17:32 2015 CET

Some references:

healpy documentation

Healpy tutorial

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:

In [2]:
for NSIDE in 2.**np.arange(6):
    print 'The number of pixels for NSIDE = %2d is %5d' %(NSIDE, hp.nside2npix(NSIDE))
The number of pixels for NSIDE =  1 is    12
The number of pixels for NSIDE =  2 is    48
The number of pixels for NSIDE =  4 is   192
The number of pixels for NSIDE =  8 is   768
The number of pixels for NSIDE = 16 is  3072
The number of pixels for NSIDE = 32 is 12288

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:

In [3]:
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'
The numbering of the pixels for NSIDE = 1 is: 
[ 0  1  2  3  4  5  6  7  8  9 10 11]

 array theta =  [ 0.26772047  0.26772047  0.26772047  0.26772047  0.5         0.5         0.5
  0.5         0.73227953  0.73227953  0.73227953  0.73227953] radians

 array phi =  [ 0.25  0.75  1.25  1.75  0.    0.5   1.    1.5   0.25  0.75  1.25  1.75] 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

In [4]:
NSIDE = 1
hp.ang2pix(NSIDE, np.pi/2, np.pi)
Out[4]:
6

Let's see what happens now if we use two separate arrays with the angular coordinates:

In [5]:
# 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)
Out[5]:
array([ 0,  1,  6, 11])
In [6]:
# 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
Atttention:  operands could not be broadcast together with shapes () (4,) (5,) 

In [7]:
# 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)
Out[7]:
array([[ 0,  1,  2,  3,  3],
       [ 0,  1,  2,  3,  3],
       [ 4,  5,  6,  7,  4],
       [ 8,  9, 10, 11, 11]])

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:

In [8]:
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":

In [9]:
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

In [10]:
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);
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.

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:

In [11]:
from astroML.datasets import fetch_wmap_temperatures
wmap_unmasked = fetch_wmap_temperatures(masked=False)
wmap_unmasked.shape
NSIDE = 512
ORDERING = NESTED in fits file
Ordering converted to RING

Out[11]:
(3145728,)

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

In [12]:
hp.nside2npix(512)
Out[12]:
3145728

Let us now proceed to its graphical representation:

In [13]:
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

In [14]:
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);
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.

No hay comentarios:

Publicar un comentario