n-d Projections !¶
hpproj allow for n-d projections, for d=3 to d=0
Caution
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.grid()
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
Out[1]:
<Table length=1>
brigthness background n_pix
MJy / sr MJy / sr
float64 float64 int64
------------- ------------ -----
2.51999285723 1.2330693081 151