`n-d` Projections !
===================

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

.. caution::

    :py:mod:`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 :

.. code-block:: python

    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

.. code-block:: python

    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 :class:`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

.. nope plot:: pyplots/proj3d.py
.. image:: pyplots/proj3d.png


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

.. code-block:: python

  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)


.. nope plot:: pyplots/proj2.5d.py
.. image:: pyplots/proj2.5d.png


2-d Projections
---------------

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

.. code-block:: python

    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 :class:`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`


.. nope plot:: pyplot/proj2d.py
.. image:: pyplots/proj2d.png


1-d Projections
---------------

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

.. code-block:: python

    from hpproj import hp_profile

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


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

.. nope plot:: pyplots/proj1d.py
.. image:: pyplots/proj1d.png


0-d Projections
---------------

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

.. code-block:: python

    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 :class:`astropy.table.Table` with the aperture photometry

.. code-block:: python

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