healpix_utils.py

#!/usr/bin/python
"""
General HealPix utils.
"""
import numpy as np
import healpy as hp

ORDER = 3
NSIDE = hp.order2nside(ORDER)
NPIX = hp.nside2npix(NSIDE)
AREA = hp.nside2pixarea(NSIDE, degrees=True) * 3600.  # in arcmin^2

def order2npix(order):
    return hp.nside2npix(hp.order2nside(order))

def get_cell(ra, dec):
    """
    Get HEALPIX cell for a given coordinates.
    """
    pos = (np.radians(dec + 90.), np.radians(ra))
    return hp.ang2pix(NSIDE, pos[0], pos[1])


def get_healpix(ra, dec):
    """
    Get position of HEALPIX pixel center.
    """
    pix = get_cell(ra, dec)
    pos = hp.pix2ang(NSIDE, pix)
    ra = np.degrees(pos[1])
    dec = np.degrees(pos[0]) - 90.
    return ra, dec


def cell_to_ra_dec(cell):
    """
    Get coordinates of the center for a given HEALPIX cell.
    """
    theta, phi = hp.pix2ang(NSIDE, cell)
    return np.degrees(phi), np.degrees(theta) - 90.


def galactic_cell_to_ra_dec(cell):
    """
    Get coordinates of the center for a given HEALPIX cell for a
    galactic grid.
    """
    from astropy.units import degree
    from astropy.coordinates import SkyCoord
    glon, glat = cell_to_ra_dec(cell)
    coord = SkyCoord(glon, glat, frame='galactic', unit=degree)
    new = coord.transform_to('icrs')
    return new.ra.deg, new.dec.deg


def change_resolution(old_order, new_order, input_pixels):
    nside_old = hp.order2nside(old_order)
    nside_new = hp.order2nside(new_order)
    center = hp.pix2ang(nside_old, input_pixels)
    if old_order > new_order:
        return hp.ang2pix(nside_new, *center)
    elif new_order > old_order:
        resol = hp.nside2resol(nside_old) * 0.01
        nside_temp = hp.order2nside(old_order + 1)
        result = []
        for d in np.array([(-1, -1), (-1, 1), (1, -1), (1, 1),
                          (-1, 0), (1, 0), (0, 1), (0, -1)]):
            points = np.array([center[0] + resol * d[0],
                               center[1] + resol * d[1]])
            result.append(hp.ang2pix(nside_temp, *points))
        result = np.unique(result)
        if new_order == old_order + 1:
            return result
        else:
            return change_resolution(old_order + 1, new_order, result)
    else:
        return input_pixels