{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# healpy tutorial\n",
    "\n",
    "See the Jupyter Notebook version of this tutorial at <https://github.com/healpy/healpy/blob/master/doc/tutorial.ipynb>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose the `inline` backend of `maptlotlib` to display the plots inside the Jupyter Notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import healpy as hp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NSIDE and ordering\n",
    "\n",
    "Maps are simply numpy arrays, where each array element refers to a location in the sky as defined by the Healpix pixelization schemes (see the [healpix website](https://healpix.jpl.nasa.gov/)).\n",
    "\n",
    "Note: Running the code below in a regular Python session will not display the maps; it's recommended to use an IPython shell or a Jupyter notebook.\n",
    "\n",
    "The resolution of the map is defined by the *NSIDE* parameter, which is generally a power of 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NSIDE = 32\n",
    "print(\n",
    "    \"Approximate resolution at NSIDE {} is {:.2} deg\".format(\n",
    "        NSIDE, hp.nside2resol(NSIDE, arcmin=True) / 60\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `healpy.pixelfunc.nside2npix` gives the number of pixels *NPIX* of the map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NPIX = hp.nside2npix(NSIDE)\n",
    "print(NPIX)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same pixels in the map can be ordered in 2 ways, either RING, where they are numbered in the array in horizontal rings starting from the North pole:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = np.arange(NPIX)\n",
    "hp.mollview(m, title=\"Mollview image RING\")\n",
    "hp.graticule()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard coordinates are the **colatitude** $\\theta$, $0$ at the North Pole, $\\pi/2$ at the equator and $\\pi$ at the South Pole and the **longitude** $\\phi$ between $0$ and $2\\pi$ eastward, in a Mollview projection, $\\phi=0$ is at the center and increases eastward toward the left of the map.\n",
    "\n",
    "We can also use vectors to represent coordinates, for example `vec` is the normalized vector that points to $\\theta=\\pi/2, \\phi=3/4\\pi$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vec = hp.ang2vec(np.pi / 2, np.pi * 3 / 4)\n",
    "print(vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can find the indices of all the pixels within $10$ degrees of that point and then change the value of the map at those indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ipix_disc = hp.query_disc(nside=32, vec=vec, radius=np.radians(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = np.arange(NPIX)\n",
    "m[ipix_disc] = m.max()\n",
    "hp.mollview(m, title=\"Mollview image RING\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can retrieve colatitude and longitude of each pixel using `pix2ang`, in this case we notice that the first 4 pixels cover the North Pole with pixel centers just ~$1.5$ degrees South of the Pole all at the same latitude. The fifth pixel is already part of another ring of pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta, phi = np.degrees(hp.pix2ang(nside=32, ipix=[0, 1, 2, 3, 4]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RING ordering is necessary for the Spherical Harmonics transforms, the other option is NESTED ordering which is very efficient for map domain operations because scaling up and down maps is achieved just multiplying and rounding pixel indices.\n",
    "See below how pixel are ordered in the NESTED scheme, notice the structure of the 12 HEALPix base pixels (NSIDE 1):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = np.arange(NPIX)\n",
    "hp.mollview(m, nest=True, title=\"Mollview image NESTED\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All `healpy` routines assume RING ordering, in fact as soon as you read a map with `read_map`, even if it was stored as NESTED, it is transformed to RING.\n",
    "However, you can work in NESTED ordering passing the `nest=True` argument to most `healpy` routines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading and writing maps to file\n",
    "\n",
    "For the following section, it is required to download larger maps by executing from the terminal the bash script ``healpy_get_wmap_maps.sh`` which should be available in your path.\n",
    "\n",
    "This will download the higher resolution WMAP data into the current directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget -c http://lambda.gsfc.nasa.gov/data/map/dr4/skymaps/7yr/raw/wmap_band_iqumap_r9_7yr_W_v4.fits;wget -c http://lambda.gsfc.nasa.gov/data/map/dr4/ancillary/masks/wmap_temperature_analysis_mask_r9_7yr_v4.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wmap_map_I = hp.read_map(\"wmap_band_iqumap_r9_7yr_W_v4.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, input maps are converted to *RING* ordering, if they are in *NESTED* ordering. You can otherwise specify `nest=True` to retrieve a     map is NESTED ordering, or `nest=None` to keep the ordering unchanged.\n",
    "\n",
    "By default, `read_map` loads the first column, for reading other columns you can specify the `field` keyword.\n",
    "\n",
    "`write_map` writes a map to disk in FITS format, if the input map is a list of 3 maps, they are written to a single file as I,Q,U polarization components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.write_map(\"my_map.fits\", wmap_map_I, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization\n",
    "\n",
    "As shown above, mollweide projection with `mollview` is the most common visualization tool for HEALPIX maps. It also supports coordinate transformation, `coord` does Galactic to ecliptic coordinate transformation, `norm='hist'` sets a histogram equalized color scale and `xsize` increases the size   of the image. `graticule` adds meridians and parallels.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(\n",
    "    wmap_map_I,\n",
    "    coord=[\"G\", \"E\"],\n",
    "    title=\"Histogram equalized Ecliptic\",\n",
    "    unit=\"mK\",\n",
    "    norm=\"hist\",\n",
    "    min=-1,\n",
    "    max=1,\n",
    ")\n",
    "hp.graticule()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`gnomview` instead provides gnomonic projection around a position specified by `rot`, for example you can plot a projection of the galactic center, `xsize` and `ysize` change the dimension of the sky patch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.gnomview(wmap_map_I, rot=[0, 0.3], title=\"GnomView\", unit=\"mK\", format=\"%.2g\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`mollzoom` is a powerful tool for interactive inspection of a map, it provides a mollweide projection where you can click to set the center of the adjacent gnomview panel.\n",
    "## Masked map, partial maps\n",
    "\n",
    "By convention, HEALPIX uses $-1.6375 * 10^{30}$ to mark invalid or unseen pixels. This is stored in healpy as the constant `UNSEEN`.\n",
    "\n",
    "All `healpy` functions automatically deal with maps with `UNSEEN` pixels, for example `mollview` marks in grey those sections of a map.\n",
    "\n",
    "There is an alternative way of dealing with UNSEEN pixel based on the numpy`MaskedArray` class, `hp.ma` loads a map as a masked array, by convention the mask is 0 where the data are masked, while numpy defines data masked when the mask is True, so it is necessary to flip the mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = hp.read_map(\"wmap_temperature_analysis_mask_r9_7yr_v4.fits\").astype(np.bool_)\n",
    "wmap_map_I_masked = hp.ma(wmap_map_I)\n",
    "wmap_map_I_masked.mask = np.logical_not(mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filling a masked array fills in the `UNSEEN` value and return a standard array that can be used by `mollview`.\n",
    "`compressed()` instead removes all the masked pixels and returns a standard array that can be used for examples by the matplotlib `hist()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(wmap_map_I_masked.filled())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(wmap_map_I_masked.compressed(), bins=1000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spherical Harmonics transforms\n",
    "\n",
    "`healpy` provides bindings to the C++ HEALPIX library for performing spherical harmonic transforms.\n",
    "`hp.anafast` computes the angular power spectrum of a map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LMAX = 1024\n",
    "cl = hp.anafast(wmap_map_I_masked.filled(), lmax=LMAX)\n",
    "ell = np.arange(len(cl))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "therefore we can plot a normalized CMB spectrum and write it to disk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(ell, ell * (ell + 1) * cl)\n",
    "plt.xlabel(\"$\\ell$\")\n",
    "plt.ylabel(\"$\\ell(\\ell+1)C_{\\ell}$\")\n",
    "plt.grid()\n",
    "hp.write_cl(\"cl.fits\", cl, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaussian beam map smoothing is provided by `hp.smoothing`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wmap_map_I_smoothed = hp.smoothing(wmap_map_I, fwhm=np.radians(1.))\n",
    "hp.mollview(wmap_map_I_smoothed, min=-1, max=1, title=\"Map smoothed 1 deg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information see the [HEALPix primer](https://healpix.jpl.nasa.gov/pdf/intro.pdf)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}