n-d Projections !

hpproj allow for n-d projections, for d=3 to d=0


healpy by default change any healpix map into the RING pixelization scheme without changing its header. Be sure to read the maps with the nest=None

First let’s setup the dataset :

import numpy as np
import matplotlib.pyplot as plt
import healpy as hp

from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.utils.data import download_file
from astropy.wcs import WCS
from astropy.table import Table

from hpproj import hp_stack

# Retrieve the Planck 857 GHz all sky map
irsa_url = 'http://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/maps/'
url = irsa_url + 'HFI_SkyMap_857_2048_R2.02_full.fits'
filename = download_file(url, cache=True)

hp_data, hp_header = hp.read_map(filename, h=True, nest=None)
hp_hdu = fits.ImageHDU(hp_data, fits.Header(hp_header))
hp_hdu.header['UNIT'] = 'MJy/sr'

# Fetch the PCCS catalog
cds_url = 'http://cdsarc.u-strasbg.fr/ftp/cats/J/A+A/594/A26/fits/'
url = cds_url + 'PCCS_857_R2.01.fits'

PCCS = Table.read(download_file(url, cache=True))
# Select a few sources
PCCS = PCCS[np.abs(PCCS['GLAT']) > 30]
PCCS = PCCS[:500]

3-d projections

hpproj allow for stacking in the healpix map

coords = SkyCoord(PCCS['RA'].data, PCCS['DEC'].data, unit="deg")
pixsize = hp.nside2resol(hp_hdu.header['NSIDE'], arcmin=True) /60 /4
hdu = hp_stack(hp_hdu, coords, pixsize=pixsize, shape_out=(128, 128))

hdu is an astropy.io.fits.ImageHDU containing the stack of all the requested positions. One can also use the keep option to retrieve all the individual maps


2.5-d projections

Using some more in-depth routine of hpproj, it is possible to place 2 sources at a relative given position on a map

from astropy.visualization import ImageNormalize, HistEqStretch

from hpproj import build_wcs_2pts
from hpproj import hp_to_wcs

coord_LMC = SkyCoord("05:23:34.60", "-69:45:22.0", unit=(u.hourangle, u.deg))
coord_SMC = SkyCoord("00h52m38s", "-72:48:01", unit=(u.hourangle, u.deg))

pair_coord = (coord_LMC, coord_SMC)

relative_pos = [0.3, 0.7]
shape_out = (512, 1024)

wcs = build_wcs_2pts(pair_coord, shape_out=shape_out, relative_pos=relative_pos)
img = hp_to_wcs(hp_hdu, wcs, shape_out)

norm = ImageNormalize(stretch=HistEqStretch(img))

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=wcs)
ax.imshow(img, origin='lower', interpolation='none', norm=norm)

ax.axhline(0.5 * img.shape[0], linestyle='--', c='k', alpha=0.5)
for pos in relative_pos:
    ax.axvline(pos * img.shape[1], linestyle='--', c='k', alpha=0.5)

2-d Projections

This is the most common projection, from an healpix map to a 2D map

from hpproj import hp_project

coord = SkyCoord(141.397513059, 40.5638050454, unit="deg", frame="galactic")
pixsize = hp.nside2resol(hp_hdu.header['NSIDE'], arcmin=True) / 60 / 4
hdu = hp_project(hp_hdu, coord, pixsize=pixsize, shape_out=(128, 128))

hdu is then an astropy.io.fits.ImageHDU containing the requested region on the sky and it’s corresponding header, which can be easily plotted with, for e.g., matplotlib


1-d Projections

The 1-d projection goes from an healpix map to a intensity profile

from hpproj import hp_profile

coord = SkyCoord(202.4865871680179, 47.181795866475426, unit="deg")
hdu = hp_profile(hp_hdu, coord)

hdu is then an astropy.io.fits.ImageHDU with the profile centered on the requested coordinates


0-d Projections

The 0-d projection goes from an healpix map to an aperture photometry, of a given position

from hpproj import hp_photometry
from astropy.coordinates import SkyCoord, Angle

coord = SkyCoord(202.4865871680179, 47.181795866475426, unit="deg")
apertures = Angle(hp.nside2resol(hp_hdu.header['NSIDE'], arcmin=True) / 60, "deg") * [7, 10, 15]
result = hp_photometry(hp_hdu, coord, apertures=apertures)

result is then an astropy.table.Table with the aperture photometry

<Table length=1>
brigthness    background   n_pix
MJy / sr      MJy / sr
float64       float64      int64
------------- ------------ -----
2.51999285723 1.2330693081   151