Add HEALPix map matplotlib visualisation

[cdeil]

I remember @woodmd recently said something along the lines of

Using HEALPix for data analysis is great, the hardest part is to visualise it.

So maybe we should try to add some things (probably a combination of helper functions and documentation) to astropy-healpix to help visualise HEALPix data?

I think for now, we should limit the discussion / additions to matplotlib, because we want to propose astropy.healpix for the Astropy core package, and they already have wcsaxes in http://astropy.readthedocs.io/en/stable/visualization/index.html.

Looking around a bit what others have:

healpy has helper functions described here and here and here. I haven't used those myself, and it's not clear if we want to look there much for guidance. It looks like their goal was to have a single function that does a ton of things via keyword arguments and it's really hard to figure out what this does. I think at the core they project to certain WCS and then call plt.imshow. If this is all they do, then probably just adding a documentation page to astropy-healpix explaining how to do this with reproject and wcsaxes would be best? (see http://reproject.readthedocs.io/en/stable/healpix.html). Or we add a similar helper function, but one that plays nicely with the existing astropy.visualisation.wcsaxes somehow?

One use case for this is the MOC.plot which @tboch added to make quickview plots like here.

At http://docs.gammapy.org/en/latest/_modules/gammapy/maps/hpxcube.html#HpxMapND.plot I see that @woodmd in addition wrote a _plot_poly method that draws the HEALPix pixels as MPL polygons. And he has the equivalent of "project to a WCS and plot" as a helper method.

In https://hips.readthedocs.io/en/latest/getting_started.html#work-with-the-result we also added a method that overplots lines to outline the HEALPix grid, calling plt.plot for the HEALPix pixel corners. This was mostly added for debugging, but possibly could be a helper function to add to astropy-healpix if others also find this useful.

Using Google to search for "plot HEALPix" another thing I found was https://github.com/pmelchior/skymapper by @pmelchior . From a quick look, I think that there are two methods in plotHealpix (see here), one calling plt.scatter (to draw a point for each pixel?) and one making a PolyCollection possibly similar to what @woodmd implemented in gammapy.maps. It's MIT-licensed, i.e. compatible with Astropy core, i.e. we can look and use if we want.

@astrofrog @woodmd @tboch @pmelchior - Thoughts on what (if any) HEALPix visualisation to add to astropy-healpix?

[pmelchior]

I made skymapper for surveys like LSST and DES, so the only currently supported map projections are conics. That's good (I'd even say: optimal) for surveys that are larger in East-West direction than in North-South. I have not coded up an all-sky projection, but I'd be willing to add it.

That said, there's a polygon and a point-per-cell method for healpix maps. The last is justified for equal-area projections (the default).

[cdeil]

@pmelchior - One thing I don't understand is why you in skymapper and the healpix guys are coding up projections at all. Why not just use astropy.wcs.WCS which is the Python wrapper for WCSlib and has all the projections already. And @astrofrog even wrote some to / from HEALPix code in http://reproject.readthedocs.io/en/stable/healpix.html a few years ago. For astropy.healpix at least we should use that for the projection-based plotting method, no?

[pmelchior]

In part, lack of information. But also, you have to make a few calculations in the projected frame to determine the form of the graticules. skymapper uses matplotlib artist shapes for that, so you need their parametric form. Simply allying the transformation isn't sufficient, you still need to understand what shape the meridians and parallels have (e.g. conics have straight meridians and circular parallels, Mollweide has straight parallels, etc.)

[woodmd]

I think it would be a great idea to have some HEALPix visualization methods built into astropy/astropy-healpix. PolyCollection is the most flexible approach but doesn't really scale well when you have a large number of pixels (> 10000). Note that the logic is also rather complex when you try to deal with all-sky maps where pixels wrap around the edges of the projection. The implementation in gammapy.maps only works for all-sky maps when the projection is centered at (0,0).

For larger maps a better approach would be a method that draws the pixels using low-level matplotlib commands. I experimented with implementing a plotting method with pcolorfast which can be used to draw sequences of arbitrary quadrilaterals but didn't succeed in getting a working example. The interface expects a sequence of adjacent quadrilaterals but for a HEALPix map there doesn't seem to be any simple way of building such a sequence. I don't know enough about the internals of matplotlib to know if this aspect of the pcolorfast interface is a hard requirement imposed by the backend. Perhaps we could propose adding a more flexibile interface to pcolorfast or we could try making our own implementation using pcolorfast as a starting point.