{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting Examples\n",
    "\n",
    "One of the core features of HPMOC is its robust and convenient suite of plotting tools.\n",
    "This page includes several examples showcasing how to use HPMOC for\n",
    "sophisticated plots. An emphasis is made on showing how to navigate different\n",
    "tools and projections, and how to quickly iterate to get a perfect plot for\n",
    "your purposes.\n",
    "\n",
    "HPMOC uses the `astropy.wcs.WCS` implementation of FITS projections for plotting. Details of the projection definitions and their associated FITS header keywords are given in [Calabretta and Greisen (2002)](https://ui.adsabs.harvard.edu/abs/2002A&A...395.1077C).\n",
    "\n",
    "## Setup\n",
    "\n",
    "We'll start by importing packages and picking some skymap files of various types from the test skymaps included in the HPMOC source repository. The loading and plotting instructions in this document work with any skymaps produced by any of the major [LIGO](https://www.ligo.org/) pipelines (in particular, anything downloaded from [GraceDB](https://gracedb.ligo.org/) should work)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from hpmoc import PartialUniqSkymap\n",
    "from hpmoc.plot import get_wcs, plot, gridplot\n",
    "\n",
    "DATA = Path(\".\").absolute().parent.parent.parent/\"tests\"/\"data\"\n",
    "\n",
    "NUNIQ_FITS = DATA/'S200105ae.fits'\n",
    "NUNIQ_FITS_GZ = DATA/'S200105ae.fits.gz'\n",
    "BAYESTAR_NEST_FITS_GZ_1024 = DATA/'S200316bj-1-Preliminary.fits.gz'\n",
    "CWB_NEST_FITS_GZ_128 = DATA/'S200114f-3-Initial.fits.gz'\n",
    "CWB_RING_FITS_GZ_128 = DATA/'S200129m-3-Initial.fits.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Plotting\n",
    "\n",
    "We can easily load and display a skymap using the default plot settings, which produce a Mollweide plot, using [hpmoc.partial.PartialUniqSkymap.plot](../hpmoc.partial.rst#hpmoc.partial.PartialUniqSkymap.plot) (a thin wrapper around [hpmoc.plot.plot](../hpmoc.plot.rst#hpmoc.plot.plot)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "m = PartialUniqSkymap.read(NUNIQ_FITS, strategy='ligo')\n",
    "# Replace the above line with the following (the GraceDB URL of the same skymap)\n",
    "# if you don't have the test files handy:\n",
    "# m = PartialUniqSkymap.read(graceid, 'bayestar.multiorder.fits,0',\n",
    "#                            strategy='gracedb')\n",
    "m.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can modify various properties; for example, we can double the default pixel height (vertical resolution) of the rendered image by doubling the `height`; the pixel `width` will not change, and so the aspect ratio will change:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(height=360)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also change the size of the figure itself, either by creating a new `matplotlib.figure.Figure` and passing it as the `fig` argument, or by simply passing the same keywords we would use for `Figure` as the `fig` argument to `plot`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(fig={'figsize': (12, 8)}, height=360, width=720)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also add color bars with the `cbar` argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(cbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot all-sky fixed HEALPix order skymaps with implicit pixel indices by passing the data arrays to `plot` as well. We can do this by instantiating a [PartialUniqSkymap](../hpmoc.partial.rst#hpmoc.partial.PartialUniqSkymap) using the pixels and ordering:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import Table\n",
    "\n",
    "t1 = Table.read(CWB_RING_FITS_GZ_128)\n",
    "PartialUniqSkymap(t1['PROB'].ravel(), t1.meta['ORDERING']).compress().plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or by passing the arguments as a tuple to [hpmoc.plot.plot](../hpmoc.plot.rst#hpmoc.plot.plot), which does the same thing behind the scenes for us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t2 = Table.read(CWB_NEST_FITS_GZ_128)\n",
    "plot((t2['PROB'].ravel(), t2.meta['ORDERING']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, this worked on both a RING-ordered and a NESTED-ordered skymap:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(t1.meta['ORDERING'], t2.meta['ORDERING'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One advantage of loading into [PartialUniqSkymap](../hpmoc.partial.rst#hpmoc.partial.PartialUniqSkymap) first is that we can call the [compress](../hpmoc.partial.rst#hpmoc.partial.PartialUniqSkymap.compress) method to eliminate redundant pixels first. This will not have a tremendous effect if we are doing a one-off plot, but it can lead to multiple orders of magnitude of performance increase (particularly with BAYESTAR skymaps) if we are going to work with the skymaps repeatedly. If you load a skymap with [PartialUniqSkymap.read](../hpmoc.partial.rst#hpmoc.partial.PartialUniqSkymap.read), as we did [earlier](#Basic-Plotting), all of this is handled for you; the skymap is even read in small chunks to avoid ever using much of a memory footprint. It also takes care of the slight differences between LIGO pipelines' skymap formats for you. Nonetheless, if you come across a HEALPix skymap that's not already supported, you can still load and plot it in HPMOC.\n",
    "\n",
    "There are many other arguments available to tweak the look and feel of the plot; see some of the specific examples below as well as the full documentation in [hpmoc.plot.plot](../hpmoc.plot.rst#hpmoc.plot.plot)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projections\n",
    "\n",
    "There are, however, many other projections available. The ones with built-in HPMOC support are discussed in this section.\n",
    "\n",
    "### Cylindrical and Pseudocylindrical projections\n",
    "\n",
    "#### Mollweide Projection\n",
    "\n",
    "The Mollweide (MOL) projection is an example of a *pseudocylindrical* projection. We can explicitly select the Mollweide projection by passing `\"Mollweide\"`, or its FITS projection code, `\"MOL\"` (or any of the aliases specified in [hpmoc.plot.get_wcs](../hpmoc.plot.rst#hpmoc.plot.get_wcs)), as the `projection` argument to `plot`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='Mollview')  # same as projection='MOL'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hammer-Aitoff Projection\n",
    "\n",
    "Other useful pseudocylindrical projections include the Hammer-Aitoff projection (AIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='AIT')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cylindrical equal-area\n",
    "\n",
    "the cylindrical equal-area projection, which preserves areas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='CEA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plate Carée Projection\n",
    "\n",
    "and the (Cylindrical) plate carée (CAR), AKA Cartesian, projection, which maps latitude and longitude coordinates one-to-one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='CAR')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zenithal Projections\n",
    "\n",
    "Zenithal projections are radially symmetric in their distortion and tend to be better suited to plotting zoomed images.\n",
    "\n",
    "#### Slant Orthographic Projection\n",
    "\n",
    "One exception is the slant orthographic (SIN) projection, which (with default settings) looks like a globe viewed from a great distance with the image plotted on its surface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='SIN')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `WCSAxes` don't handle the tick labels well for slant orthographic projections; you will get warnings and missing or strangely-positioned tick labels, as in the plots on this page, when plotting this projection. You **can** fix the labels manually by playing with `ax.coord`, but the easiest solution is to increase the plot resolution (at the expense of plotting speed), as is done in all of the remaining SIN skymaps in this document using the `height` and `width` parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(height=720, width=720, projection='SIN')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gnomonic Projection\n",
    "\n",
    "The gnomonic (TAN) projection is most often uzed for zoomed views of a specific region of the sky. Its default zoom level in HPMOC is chosen to be the same as that used by `healpy.gnomview`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='TAN')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Zenithal Equidistant Projection\n",
    "\n",
    "The zenithal equidistant (ARC) projection is similarly suitable for this purpose, though by default it produces an all-sky view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='ARC')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Zenithal Equal-Area Projection\n",
    "\n",
    "as does the zenithal equal-area (ZEA) projection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='ZEA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arbitrary `WCS` Projections\n",
    "\n",
    "These are likely to be the most useful projections for most purposes. However, any valid `astropy.wcs.WCS` instance capable of being plotted on an `astropy.visualization.wcsaxes.WCSAxes` can be passed as the projection argument. This is, in fact, how [hpmoc.plot.plot](../hpmoc.plot.rst#hpmoc.plot.plot) works behind the scenes; you can separate this step out explicitly by calling [hpmoc.plot.get_wcs](../hpmoc.plot.rst#hpmoc.plot.get_wcs) to construct a `WCS` instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = get_wcs('MOL')\n",
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which can be used as the `projection` argument. Note that [hpmoc.plot.plot](../hpmoc.plot.rst#hpmoc.plot.plot) will not automatically select the frame around the plot when passed a `WCS` instance as the projection, so we need to explicitly specify an elliptical frame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection=w, frame_class='elliptical')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Being able to use arbitrary `WCS` instances as the projection makes it easy to overplot HEALPix skymaps on the same axes as images from the many observatories that store data using FITS.\n",
    "\n",
    "Regardless of how a projection is specified, a `WCSAxes` instance will always be returned, and its `wcs` field will hold the resulting projection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "m.plot().wcs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rotations\n",
    "\n",
    "`hpmoc.plot.plot` puts a right ascension of 180 degrees at the center of the skymap (the default used by the LIGO collaboration for their skymaps). We can rotate the skymap to have a different center with the `rot` argument, which specifies the three Euler angles (rotating about the `ZXZ` reference axes) for rotating the skymap.\n",
    "\n",
    "### 2-Angle Rotations\n",
    "\n",
    "In practice, one usually simply wants to specify the coordinates at the center of the image. HPMOC handles this automatically for all built-in projections when provided only the right ascension and declination of the central point for the `rot` argument. For example, we can center a high-probability region at RA=60, dec=-15 by passing those as the `rot` argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(rot=(60, -15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, this works for any built-in projection and is usually what you want to do. The final angle corresponding to the FITS WCS `LONPOLE` keyword is calculated automatically so that the rotation is well-defined and the vertical direction at the center of the image points towards the north pole.\n",
    "\n",
    "### 3-Angle Rotations\n",
    "\n",
    "Sometimes, however, this isn't enough. You may need to specify a custom `WCS` for which HPMOC cannot automatically calculate `LONPOLE`, or you may want to apply an additional rotation to your plot. This section briefly discusses how `rot` is translated into FITS WCS keywords and how their interpretation varies accross different categories of projection.\n",
    "\n",
    "The details of how rotations are set vary by projection; in particular, the valid angle combinations are not always intuitive. One can figure out the valid angles by referencing the projection definitions in [Calabretta and Greisen (2002)](https://ui.adsabs.harvard.edu/abs/2002A&A...395.1077C), though in practice a little bit of trial an error usually suffices to figure out the correct choices. The three components of the `rot` argument tuple are used directly as the `CDELT1, CDELT2, LONPOLE` FITS header fields (explained in that document) for specifying a WCS.\n",
    "\n",
    "*HPMOC uses these coordinate definitions, non-intuitive properties and all, in the belief that cross-compatibility (with FITS in general and `astropy.wcs.WCS` in particular) outweigh the slight learning curve; that the unintuitive behaviors can be quickly learned with a little bit of documentation; and that the exceptional cases where `LONPOLE` need be specified are too rare to justify additional code complexity.*\n",
    "\n",
    "The effect of the `rot` arguments is mostly within a certain category of projections; we will therefore discuss rotations in (pseudo)cylindrical and zenithal projections with separate examples.\n",
    "\n",
    "#### Rotations in Cylindrical and Pseudocylindrical Projections\n",
    "\n",
    "*In this section we will consider the Mollweide projection, though the conclusions should mostly apply to the other* [(pseudo)cylindrical projections](#Cylindrical-and-Pseudocylindrical-projections) *(Hammer-Aitoff and plate carée) as well.*\n",
    "\n",
    "We can center the skymap on a right-ascension of 60 degrees as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(rot=(60, 0, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the Mollweide projection, we can center the plot on a point with positive declination by providing it as the second angle in `rot`. For example, we can center the plot on RA = 60, DEC = 15 as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(rot=(60, 15, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will get an error if we try to set the second Euler angle to -15 degrees, however, as we might want to do in this case to center the skymap on the region of highest probability. This can be fixed by setting the final angle to 180."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(rot=(60, -15, 180))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Rotations in Zenithal Projections\n",
    "\n",
    "*In this section we will consider the slant orthographic projection, though the conclusions should mostly apply to the other* [zenithal projections](#Zenithal-projections) *(gnomonic, zenithal equidistant, zenithal equal-area) as well.*\n",
    "\n",
    "The `rot` argument works more straightforwardly in zenithal projections, with a caveat. We can specify any right ascension and declination, e.g. a positive declination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(height=720, width=720, projection='SIN', rot=(60, 15, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or negative one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(height=720, width=720, projection='SIN', rot=(60, -15, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but, by default, the image will be rotated. Fortunately, the final `rot` tuple (corresponding to FITS `LONPOLE`) simply rotates about the central reference axis, and the image can be oriented right-side-up by setting it to 180 (the default `LONPOLE` value HPMOC uses for zenithal projections):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(height=720, width=720, projection='SIN', rot=(60, -15, 180))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Again,* [2-angle rotations](#2-Angle-Rotations) *do this automatically;* we can therefore make this same plot with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(height=720, width=720, projection='SIN', rot=(60, -15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zoomed Plots\n",
    "\n",
    "It's possible to zoom in on an area of interest by modifying the `hdelta` and `vdelta` parameters, which override the angular widths of the reference pixels. If we're zooming in on an area of interest, however, we'd be better served using a zenithal projection, e.g. the gnomonic projection (equivalent to that used by `healpy.gnomview`). We can choose a different projection from one of those specified in [`hpmoc.plot`](../hpmoc.plot.rst#hpmoc.plot.plot), specifying it by name (or one of its aliases), saving us the trouble of setting up a custom world coordinate system (`astropy.wcs.WCS`) instance. A gnomonic plot, centered at the default location with default pixel sizes (`hdelta` and `vdelta` automatically calculated), looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='TAN', rot=(60, -15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this zoom level, we can clearly see the individual HEALPix pixels, but it's hard to see the larger structure of the central probability region. We can zoom out by tweaking `vdelta` and `hdelta`, which in practice is easiest to do by trial and error until your plot looks how you'd like it to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(projection='TAN', vdelta=0.3, hdelta=0.3, rot=(60, -15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed Angular Width\n",
    "\n",
    "You can also figure out *exactly* what `hdelta` and `vdelta` mean by, again, referencing [Calabretta and Greisen (2002)](https://ui.adsabs.harvard.edu/abs/2002A&A...395.1077C) and looking up the relevant projection's interpretation of `CDELT1` and `CDELT` (the width and height, respectively, of the reference pixel in a given projection, which, it should be noted, does *not* correspond to the angular width of that pixel in most cases).\n",
    "\n",
    "This is usually not worth the trouble, but can sometimes be useful, particularly when automating plotting, though in this case the problem can sometimes most easily be solved by choosing an appropriate projection; for example, if you want a circular window of radius 5 degrees, you can easily plot it using the zenithal equidistant projection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hpmoc.plot import ARC_CDELT_BASE\n",
    "\n",
    "radius_deg = 5\n",
    "delta = ARC_CDELT_BASE * min(radius_deg, 180) / 180\n",
    "m.plot(projection='ARC', hdelta=delta, vdelta=delta,\n",
    "       rot=(60, -15), frame_class='elliptical')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scatterplots\n",
    "\n",
    "You can plot point sources, including their error regions, by passing collections of points (as [hpmoc.points.PointsTuple](../hpmoc.points.rst#hpmoc.points.PointsTuple) instances) to `plot`. You can specify each point source using its right ascension, declination (both in degrees), (optional) error region size (in degrees), and (optional) label, i.e. the name of that source (unnamed sources will be labeled with their index within the `PointsTuple`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from hpmoc.points import PointsTuple, Rgba\n",
    "\n",
    "pts1 = PointsTuple(\n",
    "    [\n",
    "        (271., 31., 2., 'One'),\n",
    "        (355., -44., 3.),\n",
    "        (15., -62., 5.),\n",
    "        (16., -60., 3.),\n",
    "        (10., -30., 10.),\n",
    "    ]\n",
    ")\n",
    "\n",
    "m.plot(pts1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The font is a little small for the point source labels; we can fix this by passing keyword arguments for `WCSAxes.text` as the `scatter_labels` argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(pts1, scatter_labels={'fontsize': 16})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also increase the DPI of the figure to increase legibility:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(pts1, fig={'dpi': 100})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can plot multiple sets of points, which can be distinguished by color, marker, and legend label (M1 location from [here](http://www.astropixels.com/messier/messiercat.html)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from astropy.units import Unit\n",
    "\n",
    "pts2 = PointsTuple(\n",
    "    points=[\n",
    "        (179., 16., 2.),\n",
    "        (62, -14, 3., 'Candidate'),\n",
    "        (\n",
    "            (5*Unit('hourangle')+34.5*Unit('arcmin')).to('deg').value,\n",
    "            (22*Unit('deg')+1*Unit('arcmin')).to('deg').value,\n",
    "            (6*Unit('arcmin')).to('deg').value,\n",
    "            'M1'\n",
    "        )\n",
    "    ],\n",
    "    rgba=Rgba(0.6, 0.2, 1, 0.5),\n",
    "    marker='2',\n",
    "    label='More points',\n",
    ")\n",
    "m.plot(pts1, pts2, fig={'dpi': 100})\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also bundle point sources with a related skymap by appending them to the [point_sources](../hpmoc.partial.rst#hpmoc.partial.PartialUniqSkymap.point_sources) field of a `PartialUniqSkymap`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.point_sources = [pts1, pts2]\n",
    "m.plot(fig={'dpi': 100})\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we can [rotate](#Rotations) and [zoom in](#Zoomed-Plots) on a region of interest; our point sources will now be plotted automatically (since they are contained in the `point_sources` field):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(fig={'dpi': 120}, rot=(70, 5), vdelta=0.3, hdelta=0.3, projection='TAN')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contours around Credible Regions\n",
    "\n",
    "We can also add contours around credible/containment regions of the plot. For example, to get the smallest region containing 90% of the GW probability, run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(fig={'dpi': 100}, cr=[.9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, with our zoomed plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(fig={'dpi': 120}, rot=(70, 5), vdelta=0.3, hdelta=0.3, projection='TAN', cr=[0.9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily add the 50% contour as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(fig={'dpi': 120}, cr=[0.5, 0.9])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.plot(fig={'dpi': 120}, rot=(70, 5), vdelta=0.3, hdelta=0.3, projection='TAN', cr=[0.5, 0.9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our example \"Candidate\" event is clearly within both contours.\n",
    "\n",
    "## Multi-layered plots\n",
    "\n",
    "Sometimes it is desirable to overlay several images. For example, we can visualize both the skymap itself and the size of its pixels by plotting them as a semi-transparent overlay using the `alpha` keyword and the same axis as the original skymap. We can include colorbars for both and pick colormaps which are visually distinguishable. We also tweak the figure size and `cbar` keyword arguments to improve the layout of the final result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = m.plot(fig={'figsize': (6, 5)}, cbar={'pad': 0.05}, cmap='bone_r')\n",
    "m.orders(as_skymap=True).plot(ax=ax, alpha=0.15, cmap='rainbow', cbar={'pad': 0.1})\n",
    "ax.set_title('S200105ae')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get rid of the point sources again, since they are visually distracting, and increase the opacity of the pixel sizes (`ORDER`) to emphasize that layer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.point_sources = []\n",
    "ax = m.plot(fig={'figsize': (6, 5)}, cbar={'pad': 0.05}, cmap='bone_r')\n",
    "m.orders(as_skymap=True).plot(ax=ax, alpha=0.4, cmap='rainbow', cbar={'pad': 0.1})\n",
    "ax.set_title('S200105ae')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we wanted to emphasize the pixel size completely, we could plot it alone and then simply plot some contours from the original skymap over it. We can entirely hide an image (but not the contours or associated point sources) by setting `cmap=None` and label and color our contours with the `cr_format` and `cr_kwargs` arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = m.orders(as_skymap=True).plot(projection='CAR', fig={'dpi': 120}, cbar=True)\n",
    "m.plot(ax=ax, cmap=None, cr=[0.5, 0.9, 0.99],\n",
    "       cr_format=lambda l, _: f\"{l*100:2.0f}% CR\",\n",
    "       cr_kwargs={'colors': ['aquamarine', 'teal', 'blue']})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can, of course, use this approach to layer skymaps over images from arbitrary FITS sources. Let's plot a random Swift-BAT high-SNR FITS image (downloaded from `astroquery`) alongside our second set of [point sources](#Scatterplots) with only the [credible region boundaries](#Contours-around-Credible-Regions) from our skymap showing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from astropy.wcs import WCS\n",
    "from astroquery.skyview import SkyView\n",
    "\n",
    "# get a Swift-BAT high-SNR FITS image from around M1 (Crab Nebula)\n",
    "hdu = SkyView.get_images(position='M1', survey='BAT SNR 150-195')[0][0]\n",
    "ax = plt.subplot(1, 1, 1, projection=WCS(hdu.header))\n",
    "ax.imshow(hdu.data, cmap='gist_heat_r')\n",
    "m.plot(pts2, ax=ax, cmap=None, cr=[0.5, 0.9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also interpolate a standard FITS image (as long as the ``astropy.wcs.WCS`` coordinate transforms are implemented in FITS) using HPMOC. This lets us put images in arbitrary projections into the same pixelization, making it easy to plot them on the arbitrary axes, do analyses with them, etc. You can see that the results are indistinguishable from the original:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interpolate the FITS image\n",
    "mh = PartialUniqSkymap(hdu.data, WCS(hdu.header))\n",
    "\n",
    "# plot a comparison\n",
    "fig = plt.figure(figsize=(8, 4))\n",
    "axw = fig.add_subplot(1, 2, 1, projection=WCS(hdu.header))\n",
    "axw.imshow(hdu.data, cmap='gist_heat_r')\n",
    "axh = mh.plot(fig=axw.figure, subplot=(1, 2, 2), projection=WCS(hdu.header))\n",
    "axh.grid(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More usefully, we can do an all-sky plot of an arbitrarily-projected FITS image with almost no effort. We can label the missing regions as gray and increase the resolution of the image to capture the small bright region in the center:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mh.plot(missing_color='gray', width=1000, height=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that by default, nearest neighbor interpolation will be used (at around 5-10 times the resolution of the original image), and units (and particularly area calculations/conversions) will *not* be handled automatically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subplots\n",
    "\n",
    "You can plot to subplots using the subplot argument. For example, we can manually make a grid of subplots using `plt.subplot` and pass the corresponding axes to [plot](../hpmoc.plot.rst#hpmoc.plot.plot). We can accomplish the same thing by passing the `subplot` argument to `plot`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2 = PartialUniqSkymap.read(BAYESTAR_NEST_FITS_GZ_1024, strategy='ligo')\n",
    "ax1 = m.plot(fig={'figsize': (16, 4)}, subplot=(1, 2, 1), cr=[0.9])\n",
    "ax1.set_title('S200105ae')\n",
    "ax2 = m2.plot(fig=ax1.figure, subplot=(1, 2, 2), cr=[0.9])\n",
    "ax2.set_title('S200316bj')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This also works with `matplotlib.gridspec.SubplotSpec` as the `subplot` argument (if you want to plot to a specific axes on a `matplotlib.gridspec.GridSpec`). We can also get the axis manually with `matplotlib.figure.Figure.add_subplot` and pass it as the `ax` argument.\n",
    "\n",
    "## Grid Plots\n",
    "\n",
    "It is often useful to compare skymaps, or look at the same skymap with different sets of point sources, plot settings, or projections. The quickest way to do this using HPMOC is to use [hpmoc.plot.gridplot](../hpmoc.plot.rst#hpmoc.plot.gridplot), which uses `GridSpec` behind the scenes. We can quickly do the same as above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.gridplot(m2, cr=[0.9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also specify multiple projections for each skymap:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.gridplot(m2, projections=['ARC', 'ZEA', 'TAN'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify how many columns (where different projections of the same skymap are grouped into a single column) appear on each row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.gridplot(m2, projections=['CAR', 'ZEA', 'TAN'], ncols=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first return value is the `GridSpec` used to arrange the subplots, while the second is a list whose rows correspond to each projection in the `projections` argument and are themselves lists of axes corresponding to each skymap. In other words, `gridplot(...)[1][p][s]` is the axes containing a plot of the s$^\\text{th}$ skymap using the p$^\\text{th}$ projection. You can use this to modify specific subplots, add titles to all plots, etc:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maps = [m, *[PartialUniqSkymap.read(f, strategy='ligo') for f in\n",
    "                        (BAYESTAR_NEST_FITS_GZ_1024, CWB_NEST_FITS_GZ_128,\n",
    "                         CWB_RING_FITS_GZ_128)]]\n",
    "gs, axs = gridplot(*maps, projections=['CEA', 'TAN'])\n",
    "for *axx, title in zip(*axs, ('S200105ae', 'S200316bj-1-Preliminary',\n",
    "                              'S200114f-3-Initial', 'S200129m-3-Initial')):\n",
    "    for ax in axx:\n",
    "        ax.set_title(title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add scatter plots of points of interest using the `scatters` argument to pass a list of [PointsTuple](../hpmoc.points.rst#hpmoc.points.PointsTuple) instances for each skymap, and we can pass separate keyword arguments to [hpmoc.plot.plot](../hpmoc.plot.rst#hpmoc.plot.plot) for each subplot (following the same shape as the returned axes). We can use this, for example, to adjust the zoom level and rotation of the gnomonic projections separately from the all-sky views and to set the `vmax` values for the color maps consistently across projections. Arguments shared across all subplots can be passed as keyword arguments to [gridplot](../hpmoc.plot.rst#hpmoc.plot.gridplot), as we will do for `vmin=0`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maps = [m, *[PartialUniqSkymap.read(f, strategy='ligo') for f in\n",
    "                        (BAYESTAR_NEST_FITS_GZ_1024, CWB_NEST_FITS_GZ_128,\n",
    "                         CWB_RING_FITS_GZ_128)]]\n",
    "gs, axs = gridplot(*maps, scatters=[[pts2]]*4, projections=['CEA', 'TAN'],\n",
    "         subplot_kwargs=[[{'vmax': s.max().value} for s in maps],\n",
    "                         [{'vdelta': 0.4, 'hdelta': 0.4,\n",
    "                           'rot': (70, 5), 'vmax': s.max().value}\n",
    "                          for s in maps]], vmin=0, ncols=1, cr=[0.9])\n",
    "for *axx, title in zip(*axs, ('S200105ae', 'S200316bj-1-Preliminary',\n",
    "                              'S200114f-3-Initial', 'S200129m-3-Initial')):\n",
    "    for ax in axx:\n",
    "        ax.set_title(title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add colorbars and layer our images as we were able to do with `plot`, though this time instead of passing `ax` to our second invocation of the plotting function, we pass the returned `GridSpec` and axes list (after tweaking some of the spacing arguments to make the layout look better). We can use this approach to make quite intricate visualizations for multiple skymaps with fairly little effort:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maps = [m, *[PartialUniqSkymap.read(f, strategy='ligo') for f in\n",
    "                        (BAYESTAR_NEST_FITS_GZ_1024, CWB_NEST_FITS_GZ_128,\n",
    "                         CWB_RING_FITS_GZ_128)]]\n",
    "gs, axs = gridplot(*maps, wspace=0.2, hspace=0.05, wshrink=0.7,\n",
    "                   scatters=[[pts2]]*4, ncols=1, top=0.95,\n",
    "                   projections=['CEA', 'TAN'], height=360,\n",
    "                   subplot_kwargs=[[{'width': 720, 'vmax': s.max().value}\n",
    "                                    for s in maps],\n",
    "                                   [{'width': 360, 'vdelta': 0.2,\n",
    "                                     'hdelta': 0.2, 'cbar': True,\n",
    "                                     'rot': (70, 5), 'vmax': s.max().value}\n",
    "                                    for s in maps]],\n",
    "                   cr=[0.5, 0.9, 0.99], vmin=0., cmap='bone_r',\n",
    "                   cr_kwargs={'colors': ['aquamarine', 'teal', 'blue']})\n",
    "for *axx, title in zip(*axs, ('S200105ae', 'S200316bj-1-Preliminary',\n",
    "                              'S200114f-3-Initial', 'S200129m-3-Initial')):\n",
    "    for ax in axx:\n",
    "        ax.set_title(title)\n",
    "orders = [s.orders(as_skymap=True) for s in maps]\n",
    "gridplot(*orders, projections=axs, fig=gs, alpha=0.3,\n",
    "         cmap='rainbow', vmin=0,\n",
    "         subplot_kwargs=[[{'cbar': {'shrink': 0.5}, 'vmax': o.max()}\n",
    "                          for o in orders], [{'cmap': None}]*4])\n",
    "plt.suptitle(\"Order vs. Probability Density Near Candidate & M1\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}