Interactive HEALPix maps with Bokeh

One of the most common plots in astrophysics is the 'all sky map', which visualises the value of some quantity across the entire sky. For instance, you can visualise the number of stars per square degree catalogued by the European Space Agency's Gaia space telescope by downloading a dataset from CDS and plotting it with the Python healpy library.

These visualisations are attractive, but really not that useful. You can compare two points on the sky and tell which has a greater density of stars, but you can't answer detailed questions like 'what are the coordinates of this interesting cluster of stars?' or 'how many stars are at this precise location?' They are also quite large files, requiring a large number of image pixels to represent the shapes of the much smaller number of HEALPix pixels (which are misaligned with the x-y pixels on your screen).

I decided to address these issues by creating an interactive version of these plots using the Python Bokeh library. Rather than plotting a 2D raster image of the HEALPix map, I directly plot the individual HEALPix pixels as vector polygons. This has several advantages:

1. Each pixel will be perfectly resolved at any zoom level.
2. We can add tooltips so that information about each pixel and its value appears when you hover over the pixel.
3. The file size is smaller because you only need to store the (x,y) vertices of each pixel polygon and its value.

I needed to solve several problems to make this work. I first calculated longitude-latitude coordinates along the four edges of each pixel using the astropy_healpix package. I then used the healpy package to rotate those coordinates into the required frame (for instance, from Equatorial to Galactic) and project them using a Mollweide projection. The Mollweide projection has a discontinuity at 90 degrees longitude and so I had to split the pixel polygons that crossed that discontinuity into two separate polygons. I used the Shapely library to then simplify these polygons, preserving topology and required that the points in the simplification were always within 0.001 of the original. Finally, I stored these polygons (and multipolygons, for those polygons split by the discontinuity) in a geopandas dataframe. 

The geopandas dataframe can be easily joined with any pandas dataframe containing the values you want to use to colour each pixel, and then plotted using Bokeh patches. For comparison, here is the same data as shown above but plotted using pixel polygons.

There are several ways that I plan to improve on this in subsequent posts!

1. Give the coordinates of the pixel in the tooltip.
2. Speed up the loading of the polygons so that I can go to greater resolution.
3. Link the map to a second plot so that you can click on a polygon and see data related to it. For instance, you could click on a pixel and see the magnitude distribution of the stars in that pixel.
4. Add the ability to switch between different maps with a slider or radio buttons.  For example, you could have the ability to show the number of stars reported by each edition of the Gaia catalogue. This would only marginally increase the size of the file because you only need one extra number per pixel for each extra map, compared to the ten or so numbers required per pixel to define the polygon.

As always, the code behind this post this blog post is available on GitHub.