Apply rotation to HEALPix map in healpy

Question

I have a HEALPix map that I have read in using healpy, however it is in galactic coordinates and I need it in celestial/equatorial coordinates. Does anybody know of a simple way to convert the map?

I have tried using healpy.Rotator to convert from (l,b) to (phi,theta) and then using healpy.ang2pix to reorder the pixels, but the map still looks strange.

It would be great if there was a function similar to Rotator that you could call like: map = AnotherRotator(map,coord=['G','C']). Anybody know of any such function??

Thanks,

Alex

Answer 1

I realise that this was asked ages ago, but I was having the same problem myself this week and found your post. I've found a couple of potential solutions, so I'll share incase someone else comes upon this and finds it useful.

Solution 1: This sort of depends on the format your data is coming in. Mine came in a (theta, phi) grid.

import numpy as np
import healpy as H
map = <your original map>
nside = <your map resolution, mine=256>
npix = H.nside2npix(nside)
pix = N.arange(npix)
t,p = H.pix2ang(nside,pix) #theta, phi

r = H.Rotator(deg=True, rot=[<THETA ROTATION>, <PHI ROTATION>])

map_rot = np.zeros(npix)

for i in pix:
    trot, prot = r(t[i],p[i])
    tpix = int(trot*180./np.pi) #my data came in a theta, phi grid -- this finds its location there
    ppix = int(prot*180./np.pi)
    map_rot[i] = map[ppix,tpix] #this being the rright way round may need double-checking

Solution 2: Haven't quite finished testing this, but just came across it AFTER doing the annoying work above...

map_rot = H.mollview(map,deg=True,rot=[<THETA>,<PHI>], return_projected_map=True)

which gives a 2D numpy array. I'm interested to know how to convert this back into a healpix map...

Answer 2

I found another potential solution, after searching off and on for a few months. I haven't tested it very much yet, so please be careful!

Saul's Solution 2, above, is the key (great suggestion!)

Basically, you combine the functionality of healpy.mollview (gnomview, cartview, and orthview work as well) with that of the reproject_to_healpix function in the reproject package (http://reproject.readthedocs.org/en/stable/).

The resulting map is suitable for my angular scales, but I can't say how accurate the transformation is compared to other methods.

-----Basic Outline----------

Step 1: Read-in the map and make the rectangular array via cartview. As Saul indicated above, this is also one way to do the rotation. If you're just doing a standard rotation/coordinate transformation, then all you need is the coord keyword. From Celestial to Galactic coordinates, set coord = ['C','G']

map_Gal = hp.cartview(map_Cel, coord=['C','G'], return_projected_map=True, xsize=desired_xsize, norm='hist',nest=False)

Step 2: Write a template all-sky FITS header (as in the example below). I wrote mine to have the same average pixel-scale as my desired HEALPix map.

Step 3: Use reproject.transform_to_healpix

reproject includes a function for mapping a "normal" array (or FITS file) into the HEALPix projection. Combine that with the ability to return the array created by healpy.mollview/cartview/orthview/gnomview, and you can rotate a HEALPix map of one coordinate system (Celestial) into another coordinate system (Galactic).

map_Gal_HP, footprint_Gal_HP = rp.reproject_to_healpix((map_Gal, target_header), coord_system_out= 'GALACTIC', nside=nside, nested=False)

It comes down, essentially, to those two commands. However you'll have to make a template header, giving the pixel scale and size corresponding to the intermediary all-sky map you want to make.

-----Full Working Example (iPython notebook format + FITS sample data)------

https://github.com/aaroncnb/healpix_coordtrans_example/tree/master

The code there should run very quickly, but that's because the maps are heavily degraded. I did the same for my NSIDE 1024 and 2048 maps, and it took about an hour.

------Before and After Images------

Answer 3

This function seems to do the trick (reasonably slow, but should be better than the for loop):

def rotate_map(hmap, rot_theta, rot_phi):
    """
    Take hmap (a healpix map array) and return another healpix map array 
    which is ordered such that it has been rotated in (theta, phi) by the 
    amounts given.
    """
    nside = hp.npix2nside(len(hmap))

    # Get theta, phi for non-rotated map
    t,p = hp.pix2ang(nside, np.arange(hp.nside2npix(nside))) #theta, phi

    # Define a rotator
    r = hp.Rotator(deg=False, rot=[rot_phi,rot_theta])

    # Get theta, phi under rotated co-ordinates
    trot, prot = r(t,p)

    # Interpolate map onto these co-ordinates
    rot_map = hp.get_interp_val(hmap, trot, prot)

    return rot_map

Using this on data from PyGSM gives the following:

hp.mollview(np.log(rotate_map(gsm.generated_map_data, 0,0)))

Upon rotation of phi:

hp.mollview(np.log(rotate_map(gsm.generated_map_data, 0,np.pi)))

Or rotating theta:

hp.mollview(np.log(rotate_map(gsm.generated_map_data, np.pi/4,0)))