Apply rotation to HEALPix map in healpy

Question

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

I tried using healpy.Rotatorto convert from (l, b) to (phi, theta) and then using healpy.ang2pixto change the order of the pixels, but the map still looks weird.

It would be great if there was a function similar to Rotatorthat you would call the following: map = AnotherRotator(map,coord=['G','C']). Does anyone know about such a feature?

Thank,

Alex

Answer 1

, , , . , , - .

1. , . (, ).

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

2: , , ...

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

2D numpy. , healpix...

Answer 2

, . , !

Saul 2, , ( !)

, healpy.mollview (gnomview, cartview orthview) reproject_to_healpix reproject (http://reproject.readthedocs.org/en/stable/).

angular, , .

----- ----------

1: cartview. , . /, coord. coord = ['C','G']

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

2: FITS ( ). , ​​ , HEALPix.

3: reproject.transform_to_healpix

reproject "" ( FITS) HEALPix. , healpy.mollview/cartview/orthview/gnomview, HEALPix () (Galactic).

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

, . , , -, .

----- ( iPython + FITS) ------

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

, , . NSIDE 1024 2048, .

------ ------

Answer 3

This function seems to do the trick (fairly 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, you get the following:

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

When rotating phi:

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

Or a twist theta:

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